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ABSTRACT 


A  numerical  technique  is  formulated,  in  a  computer  program  U2DIIF,  for  the 
solution  of  flow  over  an  airfoil  executing  an  arbitrary  unsteady  motion  in  an  inviscid 
and  incompressible  medium.  The  technique  extends  the  well  known  Panel  Methods  for 
steady  flow  into  solving  a  non-linear  unsteady  flow  problem  arising  from  the 
continuous  vortex  shedding  into  the  trailing  wake  due  to  the  unsteady  motion  of  the 
airfoil.  Numerous  case-runs  are  presented  to  verify  U2DIIF  computer  code  against 
other  theoretical  and/or  numerical  methods  as  well  as  in  cases  where  limited 
experimental  data  are  obtainable  in  literatures.  These  case-runs  include  airfoils 
undergoing  a  step  change  or  a  modified  ramp  change  of  angle-of-attack,  airfoils 
executing  harmonic  oscillation  in  pitching  and  plunging  motions  and  airfoils 
penetrating  a  sharp  edge  gust. 
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THESIS  DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may 
not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made, 
within  the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and 
logic  errors,  they  cannot  be  considered  validated.  Any  application  of  these  programs 
without  additional  verification  is  at  the  risk,  of  the  user. 
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I.  INTRODUCTION 


A.  GENERAL 

In  this  thesis,  a  numerical  method  is  formulated  and  coded  in  a  FORTRAN 
computer  program,  codename  U2DIIF  (Unsteady  2-Dimensional  Inviscid 
Incompressible  Flow),  to  solve  for  the  flow  over  an  airfoil  which  is  Executing  an 
unsteady  time-dependent  motion  in  an  inviscid,  incompressible  medium. 

B.  APPROACH 

The  basic  approach  to  this  problem  is  the  extension  of  a  very  general  and 
powerful  technique,  called  Panel  Methods,  developed  by  Hess  &  Smith  [Ref.  i]  for 
steady  potential  flow  problems,  to  include  the  unsteady  motion  of  the  airfoil  that  is 
continuously  shedding  vorticity  into  the  trailing  wake.  This  vortex  shedding  process 
creates  the  non-linearity  effects  of  the  problem  in  that  the  wake  vortices  influence  the 
flow  over  the  airfoii  which  in  turn  alters  the  vortex  shedding  as  the  airfoil  proceeds  in 
time.  It  is  this  very  non-linearity  of  unsteady  flow  that  distinguishes  itself  from  the  well 
known  steady  Panel  Methods  solution  where  the  mathematical  formulation  of  the 
problem  results  in  a  set  of  N  linear  equations  in  N  unknowns  which  are  solved  easily 
with  the  standard  Gaussian  elimination  algorithm. 

The  unsteady  flow  problem  is,  however,  deprived  of  this  relatively  easy  solution 
technique.  Instead,  an  iterative  type  of  solution  is  needed  for  this  non-linear  problem. 
The  correct  mathematical  model  must  therefore  be  formulated  to  describe  the  vortex 
shedding  process  that  provides  the  mechanism  for  the  iteration  to  proceed  towards  a 
converged  set  of  solution  in  each  time  step. 

It  is  the  objective  of  this  thesis  to  develop  a  numerical  computer  program  that 
performs  this  non-linear  potential  flow  calculation  which  proceeds  step  by  step  in  time. 
At  each  time  step,  a  complete  set  of  potential  flow  solutions,  inclusive  of  the  airfoil 
pressure  distribution,  force  and  moment  coefficients,  and  the  trailing  vortex  wake 
pattern  (strengths  and  positions  of  shed  vortices),  is  obtained. 

C.  SCOPE 

The  Panel  Method  of  Hess  &  Smith,  which  utilises  both  the  distributed  sources 
and  vorticities  as  panel  singularities,  for  steady  flow  solution  is  described  in  Chapter  II. 


Chapter  III  formulates  the  mathematical  model  for  the  unsteady  flow  problem 
and  its  solution  procedures,  highlighting  the  essential  features  in  solving  the  non-linear 
problem  of  unsteady  flow. 

Chapter  IV  describes  the  computer  program  U2DIIF,  its  essential  capabilities, 
limitations  and  the  necessary  input  set-up  for  typical  case-runs. 

The  results  of  some  of  the  case-runs  are  presented  in  Chapter  V.  They  are 
compared  with  other  theoretical  and/or  numerical  methods  as  well  as  in  cases  where 
limited  experimental  data  are  obtainable  in  the  literature.  These  case-runs  include 
airfoils  undergoing  a  step  change  or  a  modified  ramp  change  in  angle-of-attack,  airfoils 
executing  harmonic  oscillation  in  both  pitching  and  plunging  motions  and  airfoils 
penetrating  a  sharp  edge  gust. 

In  the  concluding  remarks  of  Chapter  VI,  the  future  development  and  application 
potential  of  this  numerical  method  to  other  studies  of  unsteady  2-dimensional  inviscid 
incompressible  flow  are  mentioned. 


II.  STEADY  FLOW  PROBLEM  FORMULATION 


A.  FRAME  OF  REFERENCE 

Consider  a  2-dimensional  airfoil  in  motion  with  constant  linear  velocity  -  Vqq  as 
shown  in  Figure  2.1.  Using  an  (x,y)  coordinate  system  fixed  on  the  airfoil,  where  the  x- 
axis  coincides  with  the  chord  line  originating  from  the  leading  edge  towards  the  trailing 
edge  of  the  airfoil,  the  flow  in  this  frame  of  reference  is  steady.  That  is  to  say,  the  fluid 
velocity  and  pressure  in  the  flow  field  depend  only  on  the  spatial  coordinates  (x,y)  and 
not  on  time.  The  airfoil  then  appears  to  be  submerged  in  an  onset  flow  whose  velocity 
is  V-q  and  making  an  angle  of  attack,  a,  with  the  x-axis  (see  Figure  2.1). 

B.  STEADY  FLOW  PANEL  METHODS 

1.  Definition  of  Nodes  and  Panels 

The  airfoil  surface  is  divided  into  (n)  straight-line  segments,  called  panels,  by 
(n+1)  arbitrary  chosen  points,  called  nodes,  distributed  over  the  airfoil  contour  as 
shown  in  Figure  2.2.  The  panel  numbering  sequence  starts  with  panel  1  on  the  lower 
surface  at  the  airfoil  trading  edge  and  proceeds  clockwise  around  the  airfoil  contour  so 
that  the  last  panel  (panel  n)  ends  on  the  upper  surface,  also  at  the  airfoil  trailing  edge. 

Notice  that  this  numbering  sequence  dictates  that  the  airfoil  body  always  lies 
on  the  right  hand  side  of  the  i*  panel  as  one  proceeds  from  the  Ith  node  to  the  (i+  1)* 
node.  Also  the  1st  and  the  (n+l)^  nodes  coincide  at  the  trailing  edge.  It  therefore 
facilitates,  as  shown  in  Figure  2.2,  the  common  definition  of  unit  normal  vector  n  and 
the  unit  tangential  vector  t.  for  all  panels,  i.e.,  n.  is  directed  outward  from  the  body  into 
the  flow  and  is  directed  from  the  i*  node  to  the  (i+  11th  node. 

2.  Distribution  of  Singularities 

Figure  2.2  also  indicates  that  a  uniform  source  distribution  qj  and  a  uniform 
vorticity  distribution  y  are  placed  on  the  ith  panel.  The  source  strength  q;  varies  from 
panel  to  panei  whereas  the  voracity  strength  y  remains  the  same  for  all  panels.  This 
particular  choice  of  singularity  distributions  is  one  of  the  many  types  of  singularity 
combinations  (it  happened  to  be  the  pioneering  one  though)  ever  used  in  a  wide  variety 
of  the  so  called  Panel  Methods.  The  success  of  representing  the  flow  past  an  arbitrary 
shaped  airfoil  by  surface  singularity  distributions  lies  in  the  fact  that  these  singularity 
distributions  automatically  satisfy  Laplace's  equation,  the  governing  flow  equation  for 


Figure  2.1  Frame  of  Reference  for  Steady  Flow. 
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inviscid  incompressible  flow,  and  the  boundary  condition  at  the  far  field  (301  In 
addition,  the  superposition  principle  applies  to  any  linear  homogeneous  second  order 
partial  differential  equation  such  as  Lapiace's  equation.  Therefore  one  can  buiid  up  an 
overall  complicated  How  field  by  the  combination  of  simple  flows  if  the  appropriate 
boundary  conditions  on  the  airfoil  can  be  satisfied  accurately.  For  our  case  the  overall 
flow  fleid  (represented  by  the  velocity  potential  <I>)  can  be  built  up  by  three  simpie 
flows, 


®  =  <Pco  +  <PS  +  <PV 


(eqn  2.1) 


where  (p^  is  the  potential  of  the  onset  flow, 


(9X  =*  V  ^  (x  cosa  -  y  sinai 


(eon  2.2i 


tp.  is  the  velocity  potential  of  a  source  distribution  of  strength  q(s)  per  unit  length, 


.  0(S) 

<p  =  — - In  r  ds 


(eqn  2.3) 


<p  is  the  velocity  potential  of  a  vorticity  distribution  of  strength  f(s)  per  unit  length. 


r  7(s)  o  ^ 
p  =  “  f - 0  ds 

s  3  2k 


(eqn  2.4) 


The  integrals  in  Equations  2.3  and  2.4  are  performed  along  the  surface 
contour  s  and  (r,0)  are  polar  coordinates  of  any  field  point  (x,y)  measured  from  the 
airfoil  surface  at  an  arbitrary  point  as  shown  in  Figure  2.3.  The  difficult  task  of 
evaluating  these  integrals  has  been  greatly  simplified  by  our  singularity  distributions 

postulated  :o  represent  the  flow  over  'he  airfoil;  *  hat  Is.  :nstead  of  integrating  over  the 
entire  airfoii  contour,  we  integrate  on  each  panel  along  i  straight  line  where  q.  ana  y 

are  constant,  then  sum  up  the  effects  of  all  panels.  Equation  2.1  therefore  becomes, 


Figure  2.3  Potential  Evaluation  at  a  Field  Point. 
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It  can  be  seen  from  Equation  2.5  that  <I>  is  completely  defined  if  the  (n4- 1) 
unknowns,  q.  (j=  l,2,....,n)  and  y,  can  be  calculated  using  a  numerical  technique  yet  to 
be  described.  Once  the  potential  <I>  is  solved,  the  velocity  can  be  evaluated  by  taking 
grad  <D.  At  this  point  we  introduce  a  definition  of  disturbance  potential,  (p,  as  the  sum 
of  potential  due  to  both  the  source  and  vorticity  distribution, 


<p  =  <p  +  <p 


(eqn  2.6) 


Equation  2.1  therefore  reads. 


0>  =  <Poo  +  <P 


(eqn  2.7) 


The  total  velocity  vector  is  thus, 


V  =  Vd> 

total 

=  V<Ps 


+  Vq> 


(eqn  2.8) 


The  pressure  can  be  obtained  from  Bernoulli's  Equation, 


CP  = 


P-Poc 
y*  PVoo 


1  —  1  _  (  V  total  \2 

2  '  ' 

voo 


(eqn  2.9) 


Notice  that  Figure  2.3  indicates  that  the  field  point  lies  off  the  airfoil  surface, 
however,  we  are  interested  in  field  points  that  are  on  the  airfoil  surface.  In  the  case  of 
steady  flow,  the  expressions  for  VtoUJ,  and  Cp  are  the  same  for  field  points  lying  on  or 
off  the  airfoil  surface.  It  is  nevertheless  not  the  same  in  unsteady  flow,  as  will  be  seen 
in  Chapter  III.  in  that  V  (  must  include  'he  rigid  body  motion  of  the  airfoil  when  one 
evaluates  iieid  points  on  the  airfoii  surface. 

3.  Boundary  Conditions 

The  boundary  conditions  to  be  satisfied  include  the  flow  tangency  conditions 
and  the  Kutta  Condition.  The  flow  tangency  conditions  are  satisfied  at  the  exterior  mid 
points,  called  control  points,  of  all  panels  by  taking  the  resultant  velocity  at  each 
control  point  to  have  only  (V1).  but, 


(Vn)j  «  0  ,  i-  1,2 . n 


(eqn  2.10) 


where  (Vc).  and  (Vn)j  are  the  tangential  and  normal  components  of  the  total  velocity  at 
the  control  point  of  the  panel  due  to  the  free  stream  and  the  velocities  induced  by 
the  source  and  vorticitv  distributions  on  ail  the  panels,  j  .,]=  1.2 . n). 

The  Kutta  condition  postulates  that  the  pressures  on  the  upper  and  lower 
panels  at  the  trailing  edge  be  equal  in  order  that  the  flow  leaves  the  trailing  edge 
smoothly.  By  using  Bernoulli's  equation  for  steady  potential  flow,  this  pressure 
equilibrium  condition  implies  that  the  tangential  velocities  in  the  downstream  direction 
at  the  1st  and  the  n*  panel  control  points  must  be  equal.  This  fact  is  certainly 
consistent  with  the  knowledge  that  when  steady  flow  is  established,  the  total  circulation 
over  the  airfoil  does  not  change  if  the  tangential  velocities  are  the  same  at  the  trailing 
edge  panels. 


(V1)!  =  -(Vc)n  (eqn  2.11) 

If  one  could  explicitly  express  Equations  2.10  and  2.11  in  terms  of  the  unknowns 

qj  (j  =  1,2 . ,n)  and  y,  the  task  is  then  reduced  to  solving  a  linear  system  of  (n+1) 

simultaneous  equations  for  the  (n  +  1)  unknowns. 

C.  INFLUENCE  COEFFICIENTS 

1.  The  Concept  of  Influence  Coefficients 

The  numerical  technique  employed  in  Panel  Methods  to  manipulate 
equations  2. 10  and  2.11  into  an  algebraic  system  of  linear  simultaneous  equations 
involves  the  important  concept  of  influence  coefficients.  An  influence  coefficient  is 
defined  as  the  velocity  induced  at  a  field  point  by  a  unit  strength  singularity  (be  it  a 
point  singularity  or  a  distributed  singularity)  placed  anywhere  within  the  flow  field.  In 
this  case,  it  is  the  unit  strength  singularity  distribution  on  one  panel.  Recall  that 
equations  2.10  and  2.11  simply  require  the  computation  of  the  normal  and  tangential 
.•eioctty  components  it  ill  the  panel  control  points.  The  normal  components  of 
velocities  are  essential  in  satisfying  flow  tangency  conditions  while  the  tangential 
components  of  velocities  are  necessary  for  satisfying  the  Kutta  condition  as  well  as 
computing  the  pressure  distribution.  The  procedure  is  thus  to  compute,  at  the  i^  panel 
control  point,  the  velocity  components  induced  by  the  source  and  vorticitv 
distributions  on  all  the  panels,  j  (j=  1,2 . ,n),  including  the  i*  panel  itself.  Summation 
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of  all  the  induced  velocities,  separately  for  the  normal  and  tangential  components, 
together  with  the  free  stream  velocity  components  produces  ail  the  required  (Vn).  and 
<V%  x-  1.2 . n. 

2.  Notation  for  Influence  Coefficient 

We  shall  adopt  a  consistent  set  of  notation  for  the  influence  coefficients  used 
throughout  this  documentation.  It  is  so  designated  to  permit  easy  recognition  in  that 
each  influence  coefficient  contains  all  the  associated  information  one  needs.  An 
influence  coefficient  is  denoted  with  a  superscript  and  two  subsript  as  follows: 


where  x  denotes  the  type  of  singularity  involved,  we  shall  arbitrarily  use  A,  B  and  C  for 
the  uniformly  distributed  source,  uniformly  distributed  vorticity  and  point  vortex 
respectively.  The  superscript  s  is  an  indicator  telling  which  component  the  induced 
velocity  is.  The  first  subscript  p  identifies  the  field  point  where  the  induced  velocity  is 
evaluated.  The  second  subscript  a  denotes  the  particular  singularity  contributing  to  the 
induced  velocity. 

We  thus  define,  for  the  steady  flow  problem,  the  following  influence 
coefficients  : 

•  An.j  :  normal  velocity  component  induced  at  the  1th  panel  control  point  by  unit 
strength  source  distribution  on  the  j1*1  panel. 

•  AT  :  tangential  velocity  component  induced  at  the  i ^  panel  control  point  by 
unit  strength  source  distribution  on  the  j*  panel. 

•  BTj  :  normal  velocity  component  induced  at  the  i*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  j*  panel. 

•  BT  :  tangential  velocity  component  induced  at  the  1th  panel  control  point  by 
unit  strength  vorticity  distribution  on  the  j*  panel. 

3.  Computation  of  Influence  Coefficients 

The  influence  coefficients  turn  out  to  be  related,  not  surprisingly,  to  the 
geometry  of  the  airfoil  and  the  manner  in  which  the  panels  are  formed.  Specifically,  as 
aenvcu  n  ?.ef.  21.  the  A  s  and  3'|  .niluence  coefficients.  •  aue  to  .mubrmiy  distributed 
source  or  vorticity  are  functions  of: 

•  The  natural  logarithm  of  the  ratio  of  distance  from  the  i*  panel  control  point 
(the  field  point)  to  the  (j  +  1)*  and  nodes  of  the  j1*1  panel  where  singularities 
are  distributed. 


‘C's  coefficients  will  be  needed  only  for  unsteady  flow. 

23 


•  The  angle,  in  radian,  subtended  at  the  i*  panel  control  point  (the  field  point)  by 
the  fj  —  0th  and  ;:h  nodes  of  the  ;th  panel  where  singularities  are  distributed. 

•  The  trigonometry  angles  of  the  1th  and  panels. 

Referring  to  the  geometrical  quantities  indicated  in  Figure  2.4,  the 
exoressions2  for  these  influence  coefficients  are  : 


A”j  *  [  sin(0i  - 0.)  In  Ai±J_  +  cos(e.  - Q.)  p..  J ,  for  i  *  j 

rii 


,  for  i  =  j  (eqn  2.12) 


h  _  f  cirl/fl  “  ft  \  R  —  I'rtcf  —  0  \  in  Id  1 


A ::  =  -7—  [  siniOj  -  0.)  P-.  -  cos(0.  -  0.)  In  — ^ — L  ]  ,  for  i  at  j 

‘J  ‘  J  1  J 


,  for  i  =  j  (eqn  2.13) 


3n .  = 

;i 


r:  .  _  , 


- '  :csi9.  -0.)  In  — ^ — -  —  sin(0.  —  0.)  p..  ]  .  for  :  *  i 

In  1  -  :  -  1 


,  for  1  ==  j  (eqn  2.14) 


ij  *  —  l  c°s(ei  ~  9j)  P;j  +  sin(0;  -  0.)  In  ] ,  for  i  *  j 


,  for  i  —  j  (eqn  2.15) 


where  : 


ri,i  + 1  “  V  [(xmi-xj  +  I)2  +  (ymj-yj  +  ,)2] 
r..  =  V  [(xrn  -  x;)2  +  (ym.  —  v.)2] 


I  I  1  -r  1 


ymi  ■  ‘ri(yi+yi  +  1) 


2 Actual  computation  uses  =  -  B^  and  =  A"^  to  reduce  computing 


0j  =  arctan( 


0j  =  arctanf 


Pjj  =  arctan(' 


xi  -  i  “  xi 

xm;  x-  i 


—  arctan( 


ymi  ~  . 

xmj-x/ 


D.  NUMERICAL  SOLUTION  SCHEME 
1.  Rewriting  the  Boundary  Conditions 

Using  the  concept  of  influence  coefficients,  the  flow  tangency  conditions  of 
Equation  2.10  can  be  expressed  as. 


n  n 

I  (  A";,-  q,-  1  +  Y  I  Bn.tj  +  Vx  sin(a  -  0.)  -  0  ,  i-  1,2 . n  (eqn  2.16) 

The  Kutta  condition  of  Equation  2.11,  in  terms  of  influence  coefficients,  looks 


-  1 1  A>  q.  1  -  y  £  B‘  -  Vx  cos(a  -  B,) 
j=i  i=i 

=  Z  [  A1 nj  qj  ]  +  Y  I  Bcnj  +  cos(a  -  0n)  (eqn  2. 17) 

1=1  1=1 

The  negative  signs  appearing  on  the  left-hand-side  of  Equation  2.17  are  a 
direct  consequence  of  our  definition  of  unit  tangential  vector.  In  other  words,  the 
tangential  velocities  on  the  lower  surface  panels  downstream  of  the  front  stagnation 
point  have  negative  values.  This  feature  in  fact  allows  one  to  predict  the  front 
stagnation  point  by  interpolating  the  velocity  distribution  around  the  leading  edge. 

2.  Solving  the  Strengths  of  Source  and  Vorticity  Distributions 

It  is  not  difficult  at  this  stage  to  see  that  if  we  collect  the  like  terms  in 

Equation  2.17  and  expand  Equation  2.16  for  all  i's  (i=  1,2 . ,n),  these  equations 

constitute  none  other  than  a  linear  algebraic  system  of  (n+  1)  equations  as  shown  in 
the  matrix  Equation  2.18. 
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(eqn  2.18) 
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Equation  2.18  is  a  set  of  linearly  independent  equations  which  can  be  easily 
solved  by  any  standard  linear  system  solver,  one  of  which  is  the  well  known  method  of 
Gaussian  Elimination  with  Partial  Pivoting. 

3.  Computation  of  Velocity  and  Pressure  Distribution 

Once  the  q.  (j  =  1,2,.... ,m  and  y  are  soived.  the  velocities  at  ail  the  panei 
control  points  can  be  evaluated.  Oniv  the  tangential  components  exist  since  the  normal 
components  are  already  set  to  zeroes  by  the  flow  tangency  conditions. 


V 


total  _  /v t 


00 


=  (V1).  ,  i=  1,2,. ...,n 


(eqn  2.19) 


where  : 


n  n 

(yt)i  =  £  t  A‘ij  <lj  1  +  Y  Z  Btij  +  voo  cos(a-e;)  ,  i=  1,2 . n  (eqn  2.20) 

i=l  j=l 

Substituting  Equation  2.20  into  the  Cp  equation  (Equation  2.9),  the  pressure 


coefficients  at  the  i*  panel  control  point  is 
(Co);  =  1  —  (V1).2  ,  i—  1.2 . n 


(eqn  2.21) 


4.  Computation  of  Forces  and  Moments 

The  two  dimensional  aerodynamic  coefficients  of  lift  (Cf),  drag  (Cd)  and 
pitching  moment  (Cm)  about  the  leading  edge  are  computed  by  integration  of  the 
pressure  distribution  assuming  constant  Cp  exists  in  each  panel.  The  computation  is 
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first  done  by  integrating  forces  in  the  airfoil-fixed  coordinate  system  followed  by  a 
rotation  to  the  respective  lift  and  drag  directions  along  and  perpendicular  to  the  free 
stream  (Vx)  as  follows  : 


a 


C.  -  v  fCA  (x.  :  .  -X.) 

i—  p'  1  1  1"  1  l' 

i=  1 

(eqn  2.22) 

n 

C  =  -  v  (C  ).  (v.  _  -v.) 

x  -  '  pa  -  ir  1  '  i' 

i=  1 

(eqn  2.23) 

n 

cm  =  I  (Cp)i  [(xi+1  ~  Xj)  xrr,  +  (y.+ {  -  y.)  ym.  ] 
i=  L 

(eqn  2.24) 

Cd  =  Cx  cosa  +  Cy  sina 

(eqn  2.25) 

Cp  =  C,  cosa  -  Cv  sina 

v  y  * 

(eqn  2.26) 

III.  UNSTEADY  FLOW  PROBLEM  FORMULATION 


A.  OVERVIEW  OF  UNSTEADY  FLOW  MODELING 

1.  Some  Previews 

Having  fully  understood  the  Panel  Methods  formulation  and  solution  for  the 
steady  flow  problem,  one  could  then  venture  into  the  interesting  and  complicated 
unsteady  flow  case.  In  this  Chapter,  we  shall  see  how  we  could  build  the  time- 
dependency  into  the  Panel  Methods  solution  which  has  been  proven  to  be  an  useful 
and  accurate  tool  for  steady  flow.  The  approach  in  the  unsteady  flow  problem 
formulation  will  proceed,  in  general,  in  a  manner  similar  to  Chapter  II.  However,  as  we 
go  along,  we  will  pick,  up  the  highlights  of  the  essential  differences  (also  similarity  ) 
between  the  two  problems.  Additional  flow  modeling  of  the  vortex  shedding  process 
that  greatly  influences  the  numerical  solution  technique3  will  be  discussed  in  details. 

2.  Specific  Unsteady  Flow  Model 

Recall  that  in  steady  flow,  the  problem  is  considered  solved  as  soon  as  the 

airfoil  surface  singularity  distributions  of  source  and  vorticity  q.  (j=  1,2 . ,n)  and  y  are 

determined.  These  (n+  1)  unknowns  are,  however,  time  dependent  in  unsteady  flow. 
We  therefore  introduce  a  subscript  k  as  the  time-step  counter;  that  is,  we  postulate  to 
solve  the  unsteady  flow  problem  at  successive  intervals  of  time,  starting  from  tQ  =  0.  At 

each  time-step  tk(k=l,2 . ,<*>),  we  represent  the  airfoil  by  surface  singularity 

distributions  consisting  of  source  distribution  (q.)^  (j=  1,2,... .,n)  and  vorticity 
distribution  yk.  Again  the  source  strengths  vary  from  panel  to  panel  but  the  vorticity 
strength  remains  the  same  for  all  panels. 

The  overall  circulation  Tk  at  time-step  tk  is  simply  yk  multiplied  by  the  airfoil 
perimeter,  t.  Since  the  total  circulation  in  a  potential  flow  field  must  be  preserved 
according  to  the  Helmholtz's  theorem  of  continuity  of  vorticity,  any  changes  in  the 
circulation  on  the  airfoil  surface  must  be  manifested  by  an  equal  and  opposite  change 
in  vorticity  in  the  wake.  We  cail  this  the  vortex  shedding  process  and  postulate,  as 
shown  in  Figure  3.1,  that  this  shed  vorticity  takes  place  through  a  small  straight  line 
wake  element  attached  as  an  additional  panel  to  the  trailing  edge  with  uniform  vorticity 
distribution  (yw)k-  We  shall  from  now  on  refer  to  this  panel  as  the  shed  vorticity  panel, 
The  shed  vorticity  panel  will  be  established  if  its  length  Ak  and  inclination  0k,  to  x- 

3 Referring  to  the  switch  from  a  direct  scheme  to  an  iterative  scheme. 


axis  of  the  airfoil- fixed  coordinate  system,  satisfy  the  Heimhoitz's  theorem  as  follows, 


Ak  Ww)k  +  rk  =  Fk-i 

Ak  ^w'k  =  *~k-L  ~  *~k  =  ^  r'k-l  “  '‘V 


or 


(eqn  3.1) 
(eqn  3.2) 


where  and  yk.1  are  respectively  the  total  circulation  and  vorticity  strengths  which 
are  already  determined  at  a  time-steD  tk.t  before  t. . 

Let  us  project  one  time  step  ahead  to  v  we  allow  the  shed  vorticity  panel 
to  be  detached  from  the  trailing  edge  and  get  convected  downstream  as  a  concentrated 
free  vortex,  with  circulation  Ak  (Yw)k  or  —  rk,  according  to  the  resultant  local 
velocity  occurred  at  the  center  of  the  vortex  core.  At  the  same  time  a  brand  new  shed 
vorticity  panel  is  formed  for  the  new  time  step  and  the  process  is  repeated.  Therefore 
the  shed  vorticity  panel  model  provides  exactly  the  desired  communication  mechanism 
to  carry  the  solution  from  one  time-step  to  another. 

We  now  return  back  to  the  time-step  :,  and  immediately  realise  that  as  a 
result  of  tins  continuous  vortex  shedding,  there  has  been  a  series  of  shedding  processes 
occurred  prior  to  c,  that  cummulated  in  a  string  of  concentrated  core  vortices  of 

strengths  (rk_2  —  T k_l ),  (rk_3-rk_2),  (rk.4-rk,3) . and  so  on,  forming  the 

wake  pattern  behind  the  airfoil  as  shown  in  Figure  3.1 

The  presence  of  the  shed  vorticity  panel  and  the  downstream  resultant  wake 
core  vortices  do  influence  the  upstream  flow  in  in  viscid  incompressible  flow.  In 
particular  the  shed  vorticity  panel  itself  depends  on  yk  to  determine  its  distributed 
vorticity  (yw)k,  this  in  turn  causes  changes  to  (q. )k  and  yR.  Moreover,  the  downstream 
core  vortices  that  constitute  the  wake  are  convected  under  the  influence  of  the  free 
stream  and  the  singularity  distributions  on  the  airfoil  surface  panels  including  the  shed 
vorticity  panel.  The  problem  is  thus  seen  to  be  coupled  from  this  analytical 
standpoint.  Putting  this  in  simple  mathematical  terms,  the  algebriac  system  of 
equations  ''Equation  2.13),  representing  the  flow  tangency  conditions  and  Kutta 
condition  for  steady  rlow.  are  no  longer  .inear  because  the  coefficients  u.:  ..n  the  left- 
hand-side  matrix  are  not  constants  anymore.  They  are  function  of  q^  and  y  instead. 
The  presence  of  non-linearity  is  indeed  what  drives  the  solution  scheme  into  an 
iterative  type  for  unsteady  flow. 


o.  Boundary  Conditions 

We  next  investigate  whether  our  unsteady  flow  model  is  sufficiently 
represented,  before  we  could  proceed  to  search  for  a  numerical  iterative  solution,  by 
matching  the  unknowns  with  the  available  boundary  conditions  at  time-step  tk>  Recall 
that  we  have  introduced  three  more  unknowns  •.  yw'  A,  and  0,.  in  addition  to 

\.q- }.  (j  =  1,2 ,n)  and  y, .  We  have,  however,  so  far  only  identified  an  extra  boundarv 

condition,  namely  the  Helmholtz's  theorem  (Equation  3.2)  in  conjunction  with  the  flow 
tangency  conditions  at  the  n  panel  control  points  and  the  Kutta  condition  of  pressure 
equilibrium  at  the  trailing  edge  panels.  Clearly  we  are  in  deficit  of  two  additional 
conditions  before  attempting  further  endeavour  to  solve  the  entire  system.  Basu  and 
Hancock  [Ref.  3]  suggested  the  following  assumptions  : 

•  The  shed  vorticity  panel  is  oriented  in  the  direction  of  the  locai  resultant 
velocity  at  the  panel  mid  point. 

»  The  length  of  the  shed  vorticity  panel  is  proportional  to  the  magnitude  of  the 
resultant  velocity  at  the  panei  mid  point  and  the  step  size  of  the  time-step. 

Following  these  assumptions  thus  permits  us  to  formulate  two  more  boundary 
conciuons  as  follows. 


tan0k  -  dhck 

“  (VA 


(eqn  3.3) 


\-('k-'k.l)V'[(Vw)k:  +  (Vw)k2l 


(eqn  3.4) 


where  (Uw)k  and  (Vw)k  are  the  total  velocity  components  in  x  and  y  directions  of  the 
airfoil-fixed  coordinate  system. 

The  flow  tangency  conditions  are  still, 


im-X  =  0,  i*  l,2,....,n 


(eqn  3.5) 


However,  the  Kutta  condition  must  now  include  the  rates  if  change  of 

potential  at  the  trailing  edge  panels  (unsteady  Bernoulli's  equation)  which  can  be 
related  directly  to  the  rate  of  change  of  total  circulation.  By  using  a  backward  finite 
difference  approximation  for  this  rate  of  change  of  total  circulation,  we  express  the 
Kutta  condition  as  shown  in  Equation  3.6. 


[<V),]k2  -  [<V-)„]k2 


B 


V 

tv 


K 


=  2  l  (eqn  3.6) 

B.  RIGID  BODY  MOTION  AND  FRAME  OF  REFERENCE 

Consider  a  rigid  airfoil  executing  a  time-dependent  motion,  comprising  linear 
translation  and  angular  rotation  about  the  leading  edge  in  an  inviscid  incompressible 
medium.  We  can  describe  this  arbitrary  motion  at  any  time  instant  tk  as  the  vector  sum 
of  a  mean  velocity  -Vqq,  a  time  dependent  translational  velocity  ~[U(t)  i  +  V(g  j] 
and  a  rotational  velocity  —  S!(t)  where  i  &  j  are  unit  vectors  in  the  airfoil-fixed 
coordinate  system  as  showm  in  Figure  3.2. 

If  we  continue,  as  in  steady  flow,  to  determine  the  flow  with  reference  to  the  (x.vi 
coordinate  system  fixed  on  the  airfoil,  an  observer  sitting  on  this  frame  of  reference 
sees  an  unsteady  stream  velocity,  Vsiream,  made  up  by  the  vector  sum  of  a  mean 
velocity  Vx.,  a  time  dependent  translational  velocity  rUft)  i  -  Vt)  i]  and  a  rotationai 
velocity  S2(t).  Therefore  in  this  frame  of  reference,  unlike  the  previous  ease  where  the 
airfoil  is  allowed  to  move  only  with  constant  linear  velocity,  the  flow  is  still  unsteady  in 
that  Vstream  is  time  dependent. 


Stream  =  Voo  +  KUO  i  +  V(t)  j)]  +  ft(t)  (y  i  -  x  j)  (eqn  3.7) 

We  redefine  our  disturbance  potential  to  include  the  potential  contributions  <p 
and  <pcv  from  the  shed  vorticity  panel  and  the  wake  core  vortices  respectively.  Thus  : 


w 


<p  =  q>  +  <p  +  <D  +<p 

^  ~s  ~  V  ~  w  ~  cv 


(eqn  3.8) 


We  then  write  the  total  velocity  with  respect  to  this  frame  of  reference  as. 
^  total  ^stream  ^ 


P!  Notice  that  this  total  velocity  is  obviously  NOT  the  absolute  velocity  with  respect 


direction.  Wa  have  to  make  this  distinction  ciear  because  in  our  model  on  convection 
of  core  vortices,  we  break  up  the  caiulation  into  two  steps;  we  first  convect  the  core 
vertices  using  the  resuitant  absolute  velocity  with  respect  to  an  inertial  coordinate 
system,  followed  by  determining  their  positions  with  coordinates  relative  to  the  airfoil- 
ilxed  axes. 

The  unsteady  tlow  Bernoulli's  equation  for  the  pressure  coefficients  on  the  airfoil 
surface  must  be  written  with  respect  to  the  airfoil-fixed  coordinate  system  also.  Giesing 
[Ref.  4]  showed  this  to  be  written,  in  our  notation,  as  : 


C 


P 


p~p«  =  ,  _  ,  v,^i  ,2  _  - 

«pVoo2  Vqq  Voo  Voo2  St 


(eqn  3.10) 


where  vstrearn.  and  V.  are  defined  according  to  Equations  3.7  and  3.9. 

Equations  3.8,  3.9,  and  3.10  can  be  correlated  to  their  counter-parts  in  steady 
flow,  namely  2.6,  2.S  and  2.9  respectively  with  vslream  of  Equation  3.7  replacing  the 
V  yr,  in  Equation  2.3. 

C.  TIME-DEPENDENT  INFLUENCE  COEFFICIENTS 


1.  Definition  of  Time-Dependent  Influence  Coefficients 

The  influence  coefficients,  An..,  A1--,  Bn..  and  BU,  involving  the  source  and 
vorticity  distributions  described  in  Section  C  of  Chapter  II  are  still  useful.  These  are 
indeed  time-independent  coefficients  since  they  are  functions  of  geometrical  parameters 
which  are  fixed  in  our  rigid  airfoil.  Additional  influence  coefficients  involving  the  shed 
vorticity  panel  and  the  wake  core  vortices  must  be  defined.  These  coefficients  need  to 
be  computed  in  each  time  step  since  their  positions  vary  relative  to  the  airfoil-fixed 
coordinate  system.  For  that  matter,  as  will  be  made  clear  later,  those  influence 
coefficients  involving  the  shed  vorticity  panel  must  also  be  computed  in  every  iteration 
within  each  time  step  for  the  same  reasoning. 

a.  More  A's  and  B's  Influence  Coefficients 

7  Allowing  :ne  notations  used  previously  in  steady  flow,  we  define,  with  ‘he 
use  of  the  k-suoscr.pt  to  denote  time-dependency,  additional  influence  coefficients 
involving  uniformly  distributed  singularities  of  source  and  vorticity.  They  are  the  A's 
and  B's  coefficients  : 

•  (Bnjn+i)k  :  normal  velocity  component  induced  at  the  i*  panel  control 
point  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 
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•  (Bl;  :  tangentiai  velocity  component  induced  at  the  i*  panel  control 

point  by  unit  strength  vortieity  distribution  on  the  shed  vorticity  panel  at  time 

V 

•  (Axr,  _  ,  j)k  :  x-velocity  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  source  distribution  on  the  panel  at  time  tk. 

•  UV'  _  .  :  v-veiocity  comoonent  induced  at  the  shed  vorticity  panel  mid 

1  • 

point  by  unit  strength  source  distribution  on  the  j10  panel  at  time  tk. 

•  j)k  :  x-velocity  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  vorticity  distribution  on  the  j111  panel  at  time  tk> 

•  v3-vn  _  i  :  y- velocity  component  mduced  at  the  shed  vorticity  panel  mid 

point  by  unit  strength  vorticity  distribution  on  the  panel  at  time  tk. 

•  lAxh  )k  :  x-velocity  component  induced  at  the  center  of  the  h*  core 

vortex  by  unit  strength  source  distribution  on  the  j*  panel  at  time  tk> 

•  'Ayhj)k  :  y- velocity  component  induced  at  the  center  of  the  h^  core 

vortex  bv  unit  strength  source  distribution  on  the  1th  oanel  at  time  t, . 

'  -  1  ‘  K 

•  (Bxh-)k  :  x-velocity  component  induced  at  the  center  of  the  h*  core 

vertex  by  unit  strength  vorticity  distribution  on  the  panel  at  time  tk. 

•  i'3y.nj>k  ;  y-veiocity  component  induced  at  the  center  of  the  h*  core 

vortex  by  amt  strength  vorticity  distribution  on  the  panel  at  time  tfc. 

•  iBxhn  +  i)k  :  x- velocity  component  induced  at  the  center  of  the  11th  core 
vortex  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 

V 

•  (ByhnJ-j)k  :  y- velocity  component  induced  at  the  center  of  the  h*  core 
vortex  by  unit  strength  vorticity  distribution  on  the  shed  vorticity  panel  at  time 

V 

b.  New  C's  Influence  Coefficients 

The  presence  of  discrete  core  vortices  in  the  wake  requires  the  definition  of 
new  influence  coefficients  involving  point  singularity.  They  are  the  C's  coefficients  in 
our  familiar  notations  : 

•  (Cnim)k  :  normal  velocity  component  induced  at  the  i*  panel  control 

point  by  unit  strength  mth  core  vortex  at  time  tk. 

•  (Cllrrbk  :  tangential  velocity  component  induced  at  the  ith  panel  control 

coin:  bv  unit  strength  m_n  core  vortex  it  time  . 

*  “  x 

’  C‘  _  k.  x-veioc:tv  comoonent  induced  it  the  shed  vorticitv  oanei  mid 

* »  i  .iTl  .\,  rh  ^ 

point  by  unit  strength  m  core  vortex  at  time  tk. 

•  (Cyn+1  )k  :  y-velocitv  component  induced  at  the  shed  vorticity  panel  mid 
point  by  unit  strength  core  vortex  at  time  tk. 

•  (Cxhm)k  :  x-velocity  component  induced  at  the  center  of  the  h*  core 
vortex  by  unit  strength  m’h  core  vortex  at  time  tk. 


<cx„>v 


:  y-velocity  component  induced  at  the  center  of  the  h,jl  core 


vortex  by  unit  strength  m^  core  vortex  at  time  tk. 

2.  Computation  of  Time- Dependent  Influence  Coefficients 

<B\n+  j)k  and  (Bci)n  +  are  computed  exactly  the  same  way  as  B1^  and  BE 
are  computed  using  Equations  2.14  and  2.15  with  subscript  n-i- 1  replacing  j.  Similarly, 
(Axn  +  lij)k  and  (Axhj)k  are  calculated  using  Equation  2.12  with  Q.  set  to  zero  and 
subsript  i  appropriately  replaced.  Also  (Ayn+l  j)k  and  (Ayhj)k  are  calculated  using 
Equation  2.13  with  0j  set  to  zero  and  subsript  i  appropriately  replaced.  We  do  the 
same  for  (Bxn+1.)k  and  (B'Vn  +  1  .)k  using  Equations  2.14  and  2.15  respectively  and  so 
on  for  all  the  rest  of  A's  and  B's  coefficients.  The  C's  coefficients  will  be  computed 
with  different  expressions  from  those  of  A's  and  B's  because  they  are  the  velocities 
induced  by  unit  strength  core  vortex.  It  can  be  shown  easily,  from  the  geometry  of 
Figure  3.3,  that  their  expressions  take  on  the  following  forms. 


(C"m)k  =  - 


co^a  -  (8m)k)i 

2n  (r\m\ 


(eqn  3.11) 


(cUv  -  - 


271  (rim\ 


feqn  3.12) 


where  : 


<rim>k  =  V  [(xm.  -  xm)2  +  (ymrym)2] 

xm;  =  ^(Xj  +  Xj  +  j) 

ym{  =  tt(y,  +  yi  +  1) 


xm  =  x  coordinate  of  m*  core  vortex  at  time  tk 


y  =  y  coordinate  of  mu  core  vortex  at  time 


0;  =  arctan(— — — ) 


xi  +  i_xi 


vm.  ~  v 

(9m)i,  =  arctan(- — ■ — 


xm.  -y 

i  ■'n 


By  the  same  token.  (C*  _  ,  \  and  (Cs.\  are  comDUted  bv  Equation  3.11 

while  (Cyn+  j  m)k  and  (Cyhm)k  are  computed  by  Equation  3.12  if  9,  is  set  equal  to  zero 
and  the  subscript  i  appropriately  replaced. 

D.  NUMERICAL  SOLUTION  SCHEME 
1.  The  Flow  Tangency  Conditions 

The  flow  tangency  conditions  of  Equation  3.5  can  be  written  using  the 
influence  coefficients  as  follows, 


I  (  A"ii  «j>k  ]  +  rk  I  B"  [  (V!tream)i .  a,  !k 
i-1  j-1 


+  +  Sr(C*in)(rn.!-rn)i  -  o  .  >-u . a  s.m 

m  ~  1 

where  (’vstream)i  is  evaluated  by  Equation  3.7  at  the  i115  panel  control  point 

This  equation,  though  it  seems  complex,  says  nothing  more  chan  summing  to 
zero  all  the  velocity  contributions  due  to  individual  singularity’.  Substituting  i  Ywyk  from 
Equation  3.2,  collecting  like  terms  and  rearranging  the  equation  into, 


IiA“ii(qj)kI-rt[(e/Ak)(B»u+l)k-lB«  ] 

i-1  i-l 

(  -  1  <V.„«am>|lk  •  1  -  Wl  <B\n+.>k 


-  I  I  <C“m)k  CkM  -rm)  )>  •  i-  >'2 . »  («t»  3.14) 

m=  1 

2.  The  Iterative  Solution  Procedure 

Equation  3.  Id  is  arranged  in  this  form  because  ’.ve  intend  *o  solve 
<a,)k  ‘j31 1.2 . n>  in  terms  of  yk-  Expanding  Equation  3. Id  for  ;=  1.2 . n  results  in  the 

following  matrix  equation, 


[A]{q}k  =  Yk(B}k+  {C}R 


(eqn  3.15) 


Figure  3.3  Influence  Coefficients  due  to  Point  Singularities. 
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where  [  A  ]  is  an  n  *  n  matrix  whose  elements  are  known  constants.  {  B  }k  and  {  C  }k 
are  n  x  1  column  vectors  whose  elements  are  known  only  if  the  shed  vorticity  panel  at 
tk  is  established;  that  is.  if  Ak  and  0;<  are  known,  then  we  can  calculate  all  the 
infiuence  coefficients  on  the  right-hand-side  of  Equation  3.14.  We  therefore  make  use 
of  this  idea  to  formulate  our  iterative  solution  procedure  as  follows  : 

(1)  Project  the  wake  core  vortices  downstream  according  to  the  time  step  and  the 
local  resultant  velocities  at  their  respective  centers  with  respect  to  an  inertial 
coordinate  system. 

(2)  Compute  the  coordinates  of  these  core  vortices  relative  to  the  airfoil-fixed 
coordinate  system  due  to  its  time-dependent  motion. 

(3)  Stan  iteration  cycle  for  current  time  step  by  initially  assuming  some  guess 
values  of  Ak  and  0k-  We  can  use,  except  for  the  first  time-step,  values 
obtained  at  previous  time  step.  Compute  then  the  influence  coefficients  needed 
in  Equation  3.14  or  3.15. 

(4)  Obtain  (q.ik  in  terms  of  yk  by  solving  Equation  3.15  as  a  linear  system  with 
two  right-hand-stdes  by  the  same  Gaussian  elimination  algorithm  used  in 
steady  flow. 

>5)  Calculate  the  tangential  velocities  at  the  trailing  edge  panels,  all  in  terms  of 

■V 

(6)  Invoke  the  Kutta  condition  of  Equation  3.6  (with  some  efforts  in  algebriac 
manipulation)  to  solve  for  yk  since  it  is  the  only  unknown  in  that  equation. 

(7)  Once  yk  is  solved,  (qpk  are  then  known.  We  can  then  calculate  the  velocity 
components  (L'w)k  and  (Vw)k  at  the  mid  point  of  the  shed  vorticity  panel. 

(8)  Equation  3.3  and  3.4  hence  enable  us  to  update  the  values  of  Ak  and  0k. 

(9)  Repeat  the  iteration  cycle  from  steps  (3)  to  (9)  until  converged  values  of  Ak 
and  ©k  are  obtained.  Alternatively  convergence  can  be  set  for  (Uw)k  and 
(Vw)k  instead. 

(10)  Compute  the  tangential  velocities  and  disturbance  potential  at  all  panel 
control  points  in  order  to  determine  the  pressure  distribution  which  can  be 
integrated  to  give  forces  and  moments. 

(11)  Compute  the  resultant  velocities  which  occur  at  the  centers  of  all  the  core 
vortices  that  will  be  convected  down-stream.  These  velocities  must  be  the 
absolute  velocities  with  respect  to  an  inertial  coordinate  system. 

3.  Computation  of  Velocities 

The  iterative  procedure  mentioned  in  the  previous  subsection  requires 
calculation  of  tangential  velocities  at  the  trailing  edge  panels  and  the  absolute  velocity 
components  (Uw)k  and  (Vw)k.  They  are  computed  differently  due  to  the  use  of  a 
different  frame  of  reference. 


a.  Tangential  Velocities  on  Airfoil  Panels 

The  tangential  velocities  [( Vt).]k  (i=  1,2 . n)  at  all  the  panel  control  points 

are  calculated  using  the  airfoil-fixed  coordinate  frame  of  reference  as  follows  : 


KVU  -  £  I  A',  (1\  1  +  Yk  I  S', 

i-1  i=l 


+  ^  (Vsireamt  •  *i  k+  (Bl|4i+  l\ 


+  £KC'im)k(rm-l-rm)l  •  ^ . - 

m  —  1 


(eqn  3.16) 


b.  Core  Vortices  Convection  Velocities 

The  resultant  velocities  at  all  core  vortices  are  calculated  using  the  inertial 
frame  of  reference  fixed  with  respect  to  but  resolving  them  into  components  in  the 
directions  coincide  with  the  airfoil-fixed  coordinate  svstem  as  shown  below  : 


w„>k  -  1 1  (Ayk  (qj>k  i  +  ykl  (B*hiv 

i=i  i=i 

+  (voo  ■  i  )k+  (Yw)k  j)k 


+  £[(cta)k(rm-.-rm)i 


(eqn  3.17) 


(v„)k  *  I  [  <Ayk  (qpit  1  +  rkI  (uyk 

1-1  i-1 

+  <V*,.j)k+<Y„.)k<BV„+1\ 


+  £i(cvk(rm.,-rra)i 

m  =  1 
m*  h 


(eqn  3.18) 


Notice  the  use  of  instead  °fvstream  in  Equations  3.17  and  3.18.  Also 
the  subscript  h  is  just  an  index  usable  for  any  core  vortex.  We  can  obtain  (L'w)k  and 
(Vw)k  if  h  is  replaced  by  n+  1  in  these  equations. 

4.  Disturbance  Potential  and  Pressure  Distribution 
a.  Why  We  Seed  the  Disturbance  Potential 

The  concept  of  disturbance  potential  <p  has  been  instrumental  in  the 
formulation  of  both  the  steady  and  unsteady  flow  problems.  However,  it  has  never 
gone  beyond  using  it  merely  as  a  vehicle  to  understanding  the  superposition  of  simple 
flows.  The  disturbance  potential  need  not  be  solved  for  at  all  in  the  steady  flow 
problem  formulation.  This  is  because  what  one  really  is  going  after  is  the  spatiai 
derivative  of  this  disturbance  potential,  i.e.  the  disturbance  induced  velocity,  from 
which  the  pressure  distribution  can  be  obtained.  We  have,  in  all  our  solutions  so  far. 
been  successful  in  avoiding  any  disturbance  potential  calculation  since  the  concept  of 
influence  coefficients  allows  us  a  direct  evaluation  of  the  velocity.  Unfortunately,  as 
can  be  seen  in  Equation  3.10,  when  we  proceed  further  to  compute  the  pressure 
distribution  on  the  airfoil  surface  in  unsteady  flow,  we  are  faced  with  the  problem  of 
evaluating  the  disturbance  potential  <p.  or  more  precisely  the  rate  of  change  of©,  which 
we  approximate  by  using  a  backward  finite  difference  expression.  Therefore,  the 
pressure  coefficients  at  the  i*  panel  control  point  can  be  rewritten,  in  terms  of  non- 
dimensional  variables,  as, 

=  KV^lk2  -  [(VU2  -  2 (eqn  3.19) 

Tc  Tc-l 

where  (V1^  is  calculated  by  Equation  3.16  and  (Vslreain)-  is  the  non-dimensional  (by 
Vqq)  form  of  Equation  3.7  evaluated  at  the  i*  panel  control  point. 

We  thus  need  to  calculate  at  each  time  step,  the  disturbance  potential  at  all 
the  panel  control  points.  Short  of  having  to  solve  the  Laplace's  equation  by  a  finite 
inference  scheme,  we  evaluate  'he  disturbance  potential  ©  by  integrating  "he  velocity 
field  in  two  stages  from  upstream  at  infinity  to  the  airfoii  leading  cage,  then  along  the 
airfoil  surface  from  the  leading  edge  to  each  panel  control  point.  Care  must  be  taken 
here  to  include  only  the  velocity  contribution  due  to  disturbances. 

One  important  question  arises,  in  this  approach,  as  to  what  value  of 
disturbance  potential  we  should  use  at  infinity  before  we  carry  out  the  line  integral.  We 


muse  therefore  analyse  the  behaviour  of  <p  at  infinity  by  examining  the  singularities  that 
constitute  the  disturbance.  They  are  the  source  and  vorticity  distributions  on  the  airfoil 
surface  and  the  core  vortices  in  the  wake.  These  singularities  induce  no  velocity  at 
infinity  from  the  knowledge  of  simple  flows.  In  other  words,  the  disturbance  potential 
(p  at  infinity  is  independent  of  spatiai  coordinates.  The  next  question  we  should  ask  is 
whether  <p  at  infinity  is  timq-dependent?  Let  us  adopt  the  view-point  that  if  we  are  at 
infinity  looking  ac  our  airfoil  and  its  associated  wake,  we  simply  see  a  point  vortex  with 
a  total  circulation  T0  at  time  tQ.  We  have  already  identified  that  T0  remains  constant 
by  Helmholtz's  theorem.  It  only  gets  redistnbuted.  as  time  progresses,  over  the  airfoil 
surface  and  in  the  wake.  Notice  that  the  previous  statement  regarding  what  one  would 
see  at  infinity  said  nothing  about  the  source  distributions.  The  source  distributions 
though  vary  (or  get  redistributed  )  as  the  time  progresses,  the  total  source  strength 
necessarily  remains  zero  at  all  time  in  order  to  enforce  a  closed  contour  representing 
the  airfoil  thickness.  This  is  also  the  reason  why  the  unsteady  flow  solution  needs  an 
additional  model  to  handle  the  vorticity  conservation  since  the  source  conservation  is 
already  implicitly  so  for  a  closed  contour  to  exist.  From  these  discussions,  we  are 
certain  that  the  disturbance  potential  to  at  infinity  is  an  absolute  constant  (independent 
of  time  and  spatial  coordinates)  whose  value  is  fixed  only  by  the  initial  condition  we 
decide  to  start  solving  the  unsteady  problem.  The  actual  value  of  cp  at  infinity  is  in  fact 
immaterial  so  long  as  we  know  it  is  constant  because  its  value  disappears  conveniently 
as  we  subtract  (q>.)k_j  from  (<p;)k  in  Equation  3.19. 
b.  Computation  of  Disturbance  Potential 

We  begin  by  choosing  an  arbitrary  straight  line  extending  upstream  to 
infinity  from  the  leading  edge  of  the  airfoil  along  a  direction  parallel  to  Vqo.  For 
practical  purposes,  we  set  infinity  at  say  ten  chord  lengths  away  from  the  leading  edge 
since  the  velocities  induced,  at  field  points  thereafter,  by  the  disturbances  are  small 
enough  to  be  negligible.  This  line  is  divided  into  z  panels  with  element  lengths  near  the 
leading  edge  comparable  to  the  panel  sizes  used  on  the  airfoil.  However,  the  panel  size 
is  progressively  increased  to  take  advantage  of  the  inversely  decaying  induced  velocities 
at  larger  distances.  We  compute  the  tangenuai  components  of  the  induced  velocities  at 
the  mid  points  of  these  panels  using  influence  coefficients  analogous  to  those  used  on 
the  airfoil  panels.  Using  subscript  f  to  denote  these  panel  mid-points,  we  can  define 
influence  coefficients  (A‘rj)k,  (Atfjn+1)k,  (B^,  (Bcf)n  +  1)k,  and  (Ccfm)k  and  compute 
them  using  the  same  expressions  for  calculating  the  A's,  B's  and  C's  coefficients  used 


before  with  cos0.  replaced  by  ( —  cosa),  sin0;  replaced  by  ( —  sina)  and  subscript  i 
replaced  by  f.  With  the  help  of  these  influence  coefficients,  the  tangential  velocity 
induced  by  disturbances  at  the  1th  panel  mid  point  is : 


^  VfJk  Z  f  +  Ac  —  l  BVk 

1=1  1=1 


k-1 


ca  —  1 


(eqn  3.20) 


valid  for  f=  1,2,. ...,z.  The  disturbance  potential  at  the  airfoil  leading  edge  is  the  sum  of 
the  products  of  the  disturbance  induced  velocity  at  each  panel  and  the  panel  length. 


«PiA  = 


z  v  [  -xf~  i  v*  (t" f-r  i  Vf)~ 

f=  i 


(eqn  3.21) 


Similariy,  for  the  line  mtegral  over  the  airfoil  surface,  we  compute  the 
tangential  component  of  the  disturbance  induced  velocity  at  the  ith  panel  control  point 
using  the  following  equation  : 


KVV*Jk  Z  t  A\j  fapk  1  +  ^k  S  fitij 

j- I  j= 1 


k-1 


+  (YA  (BW i)k  +  I  [  -  rm)  1 

m  =  1 


(eqn  3.22) 


which  is  valid  for  i=  l,2,....,n.  Performing  the  line  integration  by  summation,  the 
disturbance  potential  at  the  i*  nodal  point  on  the  airfoil  is : 


<0 


notie  rk  -^le'k  —  ^  <ptyk  ' j , j  —  i  ’  'or  a  >  ‘  “  ‘le 

j  =  ile 

V1 


=  tyA  “  Z  [(vt<p>jlk  rj,j  +  l  *  fori!e  >  1  -  1  (eqn  3.23) 
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where  r. .  + 1  denotes  the  panel  length. 
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Finally,  the  disturbance  potential  at  the  1th  panel  control  point  is, 

(<Pi)k  =  K  [  (<Pnode  j)k  +  (<Pnode  i+  l\  5  *  i=  l’2 . ’n  (ecin  124) 

5.  Computation  of  Forces  and  Moments 

The  Cj,  Cd  and  Cm  about  the  leading  edge  are  calculated  in  exactly  the  same 
way  as  it  is  done  for  the  steady  flow  problem  by  integrating  the  pressure  distribution 
(See  section  D-4  of  Chapter  II). 

E.  FLOW  MODELING  OF  SHARP  EDGE  GUST  FIELD 

The  unsteady  flow  solution  described  so  far  can  be  extended  to  the  study  of 
airfoils  penetrating  a  sharp  edge  gust  by  modifying  the  boundary  conditions  with  the 
assumption  that  the  gust  front  remains  straight  while  passing  through. the  airfoil.  The 
same  assumption  has  been  used  in  both  [Ref.  3]  and  [Ref.  4],  An  additional  model  m 
[Ref.  3]  using  distribution  of  singularities  along  the  gust  front  had  successfully 
attempted  to  simulate  the  distortion  of  the  gust  front  passing  over  the  airfoil  surface. 
It  was  shown  that  the  pressure  distributions,  during  the  time  when  the  gust  front 
remained  on  the  airfoil  surface,  were  affected  only  at  the  neighbourhood  of  the  gust 
front.  The  overall  pressure  upstream  and  downstream  of  the  gust  front  stayed 
essentially  the  same.  The  distorted  gust  front  model  is  not  used  in  program  U2DIIF. 
The  use  of  the  relatively  simple  yet  sufficiently  accurate  model  of  a  straight  gust  front 
affords  the  modifications  to  the  unsteady  flow  solution  to  be  confined  only  to  the  flow 
tangency  conditions.  That  is  to  say,  the  expression  of  Vstream  in  Equation  3.7  would 
include  the  gust  velocity  for  panels  that  are  already  in  the  gust  field  during  the 
penetration  phase.  Similarly,  the  computation  of  core  vortex  velocities  using 
Equations  3.1"  and  3.18  have  the  gust  velocity  iuded  to  ‘he  Vx  if  the  core  vortices  ire 
already  m  the  gust  Held. 

In  an  attempt  to  generalise  the  solution  for  cases  where  airfoils  enter  the  gust 
field  at  an  angle  of  attack,  the  convenient  model  used  in  [Ref.  3]  by  setting  the 
computation  to  proceed,  for  the  undistorted  gust  front  simulation,  so  that  the  gust 
front  always  coincides  with  the  nodal  points  is  difficult  to  implement.  At  any  one  time, 


if  an  airfoil  enters  a  gust  field  at  an  angle  of  attack,  the  gust  front  would  appear  in 
between  two  nodes  of  a  particular  panel  on  one  surface  while  the  gust  front  proceeds 
from  node  to  node  on  the  other  surface.  We  therefore  further  modify  the  How 
tangency  condition  only  on  that  particular  panel  where  the  gust  front  lies  in  between 
two  nodes  by  taking  the  gust  velocity  on  that  panel  to  be  proportional  to  the  fraction 
of  panel  length  partially  submerged  in  the  gust  field. 


IV.  DESCRIPTION  OF  COMPUTER  CODE  L2DIIF 


A.  PROGRAM  U2DIIF  STRUCTURES  AND  CAPABILITIES 
1.  Restrictions  and  Limitations 

The  numerical  formulations  of  both  the  steady  and  unsteady  flow  problems 
outlined  in  the  previous  Chapters  are  coded  in  a  FORTRAN  computer  program  called 
U2DIIF  (See  Appendix  A  for  the  program  listings).  The  present  solution  methods 
treat  the  inviscid  and  incompressible  flow  as  an  approximation  :o  the  real  flow  so  long 
as  the  viscous  effect  is  negligible  and  the  flow  stays  attached  on  the  airfoil  surface  at  all 
time.  These  restrictions  are  no  strangers  to  my  me  who  :s  familiar  with  any  other 
potential  flow  sciutton  methods.  Other  'ha n  'he  .mndcit  restrictions  of  potential  flow 
solution,  the  method  is  entirely  general  in  that  the  shape  of  airfoil  is  arbitrary  and  any 
arbitrary  continuous  motion  of  the  airfoil  could  be  simulated  using  either  the  closed 
form  i.e.  explicit  equations-  or  iiscrete  data  mints  '0  eescrihe  me  'ime-histcry  ef  the 
translational  and  rotational  velocities. 

The  storage  of  the  computer  that  carries  out  the  calculations  may  be  the  other 
limitation  one  should  consider.  The  storage  requirements  grow  rapidly  with  the  number 
of  panels  (n)  and  the  number  of  computation  time  steps  um).  By  far  the  prime 
contributor  to  this  storage  requirement  comes  from  the  massive  amount  of  influence 
coefficients.  The  number  of  influence  coefficients  increases  with  the  square  of  the 
number  of  panels  (n2).  Each  time  step  increment  adds  (2n  +  m2)  more  influence 
coefficients  due  to  the  formation  of  shed  vorticity.  The  current  program  fixes  the 
maximum  number  of  airfoil  panels  to  200  and  the  maximum  allowable  time  steps  is 
also  200. 

An  additional  constraint  worth  mentioning  concerns  the  gust  field  simulation 
whereby  the  current  solution  methodology  is  valid  except  in  the  use  of  the  same 
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irrotationality  of  the  flow  field.  There  is  no  doubt  that  the  flow  fields  upstream  and 
downstream  of  the  gust  front  are  irrotational.  However,  when  one  needs  to  obtain  the 
pressure  on  the  airfoil  surface,  an  implicit  integration  is  done  across  the  gust  front.  A 
flow  field  inclusive  of  the  gust  front  is  rotational  since  the  line  integral  of  velocity  in  a 


closed  path  does  not  vanish  when  the  gust  front  is  crossed.  Failing  the  proper 
derivation  of  a  new  pressure  equation  applicable  to  unsteady  rotational  flows,  care 
must  be  exercised  to  regard  the  present  method  as  an  approximate  solution  to  gust 
fields  of  weak  strengths  only. 

2.  Current  Structures  of  U2DIIF  MAIN  Program 

The  overall  program  logic-flow  chart  is  as  shown  in  Figure  4.1.  The  program 
first  reads  in  the  input  data  from  filecode  1  and  sets  up  the  airfoil  panel  nodes  and 
slopes.  Immediately  after  that,  the  steady  flow  calculations  are  executed  for  the  initial 
angle  of  attack  a.  according  to  the  solution  scheme  described  in  Section  D  of 
Chapter  II.  The  steady  flew  solution  is  included  primarily  to  : 

•  Provide  the  necessary  initial  parameters  for  the  unsteady  flow’  solution  to 
proceed  in  time.  In  other  w'ords,  the  steady  flow  solution  handles  the  and 
initial  angle  of  attack  a.  one  decides  to  begin  the  unsteady  flow  calculation. 

•  Allow  the  code  to  function  directly  as  a  steady  flow  soiver  as  and  when 
necessary  without  having  to  do  the  time  consuming  unsteady  flow  iterative 
solution  and  approach  the  steady  flow’  as  time  approaches  infinity. 

The  program  ’ermir.ates  once  the  steady  flow  calculations  are  done  if  the 
program  determines,  based  on  the  input  data  set  oy  user,  no  requirement  for  unsteady 
flow  solutions.  Otherwise  the  unsteady  flow  calculations  wall  be  activated  by  selecting 
and  computing  the  rigid  body  motions  of  the  airfoil  and  the  corresponding 
computation  time-step  size.  Currently,  all  the  time  dependent  motions  are  equation¬ 
generated,  they  are  the  positions  and  rates  of  the  translational  and  rotational  motions. 
Incorporated  as  case-runs  within  the  program  U2DIIF  are  the  following  motions  : 

(1)  Step  change  in  angle  of  attack  from  any  initial  value. 

(2)  Modified- ramp  change  in  angle  of  attack  about  any  pivot  point  from  any 
initial  value. 

(3)  Harmonic  translational  motion  at  any  angle  of  attack. 

(4)  Harmonic  rotational  motion  about  any  pivot  point  at  any  mean  angle  of 
attack. 

(5)  Sharp  edge  gust  penetration  at  any  angle  of  attack. 
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function  of  time,  the  program  couid  be  easily  modified. 

The  computation  time-step  sizes  for  the  harmonic  translational  and  rotational 
motions  are  constant  values  determined  by  the  frequencies  (FREQ)  and  the  number  of 
computation  per  cycle  (DTS).  For  the  case  of  step  change  in  angle  of  attack,  the 
computation  time-step  size  is  progressively  increasing,  from  a  starting  value  (DTS),  as 
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Figure  4.1  Flow  Chart  for  U2DIIF  Computer  Code. 


<?.  -’I  -t.-' 


-Y"-.  C.a'.i  ‘  A..Y.1 


Time  SteD  increment 


Compute  Rigid  Body  Motions  and 
Convect  Core  Vortices  Downstream 


Convergence  of 

iVk  *  (Vic  ? 


te  Influence  Coefficients 


Solve  (qj^  in  terms  of  yk  by 
Gaussian  Elimination  with  2  R.H.3. 


Invoke  Xutta  Condition,  Solve  y  1 


Update  (U  ) 


A.  and  0. 


Compute  Velocity  [  < Vc) . J ,  Disturbance 
Potential  (<P|)k  and  Pressure  [(CU 


Compute  Aero  Coeff  (Cg)k,  (Cj).  and  (C  )k 


■ — Adjust''''-- 
Time-Step  ? 


omoute  ,-cesuxtant  eioc 


Re-initialising  Parameters  Before 
Starting  the  Next  Time  Step 


Figure  4.1  .  (cont  d.) 


v  v  W 


time  increases.  The  case  of  the  modified  ramp  change  in  angle  of  attack,  adopts 
initially  a  constant  computation  time-step  size  (DTS)  during  the  transient  rising  of  the 
angie  of  attack.  Once  the  final  angle  of  attack  is  reached,  the  computation  time-step 
size  is  progressively  increased  also.  Similarly  for  the  case  of  airfoil  penetrating  a  gust 
tie  id.  the  oomoutation  time-step  size  (DTS')  ;.s  constant  during  the  period  when  the  gust 
front  remains  on  the  airfoil  surface  but  progressively  increases  once  the  entire  airfoil  is 
submerged  in  the  gust  field.  These  variations  in  computation  time  steps  are  to  provide 
greater  flexibility  both  in  capturing  transients  and  covering  relatively  large  total  time  of 
computation  without  having  to  contend  with  the  storage  space  requirements  described 
previously.  These  variations  in  time-step  sizes  described  so  far  are  associated  with 
setting  the  input  parameter  TADJ  to  zero.  If  TADJ  is  chosen  to  be  non-zero,  all  the 
case-runs  would  compute  initially  using  the  starting  time-step  sizes,  based  on  DTS  for 
non-osciilatory  motions  and  FREQ  Sc  DTS  for  harmonic  motions,  and  the  program 
would  prompt  for  an  user  choice  of  time  step  adjustment.  If  the  answer  is  yes,  the 
program  would  back-track  the  previous  solution  and  recompute  the  current  solution 
using  an  adjusted  time-step  size  that  is  TADJ  times  the  initial  value  <DTSV  This  special 
time  step  variation  feature  gives  the  program  added  capability  of  allowing  an 
interactive  time  step  selection  during  the  progress  of  unsteady  How  computation.  The 
ability  to  back-track  and  recompute  the  current  solution  using  a  different  time-step  size 
enhances  the  possibility  of  using  program  U2DIIF  together  with  a  viscous  flow  solver 
forming  an  Inviscid- Viscous-Interactive  solution  scheme  which  often  requires  such  time 
step  variations. 

The  MAIN  program  performs  the  iterative  solution  procedures  set  out  in 
Section  D  of  Chapter  III.  The  convergence  check  during  the  iterative  solution  is  done 
through  the  user  specified  tolerance  between  successive  iterative  solutions  of  both 
(Uw)k  and  (Vw)k.  The  solution  continues  into  the  next  time-step  by  selecting  the  time 
step  size  according  to  the  particular  case-run  and  projecting  all  the  wake  core  vortices 
downstream  so  that  their  new  positions  relative  to  the  airfoil  at  the  new  time  step  can 
be  correct; v  determined. 

B.  DESCRIPTION  OF  SUBROUTINES 
1.  Subroutine  BODY 

This  subroutine  is  called  by  subroutine  SETUP  if  the  user  selects  an  airfoil 
that  is  either  a  NACA  XXXX  or  230XX  type.  It  in  turn  calls  subroutine  NACA45  to 
obtain  the  airfoil  thickness  and  camber  distributions  and  returns  with  the  computed 
(x,y)  coordinates  of  the  panel  nodal  points. 


2.  Subroutine  COEF 

This  subroutine  is  called  by  the  MAIN  program  in  the  unsteady  flow 
calculations.  It  utilises,  at  each  iteration  cycle,  the  influence  coefficients  generated  by 
subroutine  INFL  to  calculate  the  coefficients  of  the  matrix  Equation  3.15  by  expanding 
Equation  3.14.  These  matrix  coefficients  are  necessarily  set  up  in  this  way  so  that  the 
source  strengths  could  be  solved  'in  terms  of  the  vorticity  strength  by  subroutine 
GAUSS  as  a  linear  system  with  two  right-hand- sides. 

3.  Subroutine  CO  FISH 

This  subroutine  is  called  by  the  MAIN  program  to  set  up  the  coefficients  of 
the  matrix  system  of  Equation  2.18  for  steady  flow  where  the  source  strengths  and 
vorticity  strength  are  solved  simultaneously  by  subroutine  GAUSS  as  a  linear  system 
with  one  right-hand-side.  The  matrix  coefficients  are  calculated  using  Equations  2.16 
and  2.i". 

4.  Subroutine  CORVOR 

This  subroutine  is  called  by  the  MAIN  program  at  nearing  the  end  of  the 
unsteady  riow  calculations  before  starting  a  new  time  step.  It  computes  the  resultant 
convective  velocities  for  all  the  wake  core  vortices  with  respect  to  an  inertial  frame  of 
reference  according  to  Equations  3.1"  and  3.18  where  all  the  appropriate  influence 
coefficients  are  linked  through  common  block  from  subroutine  INFL. 

5.  Subroutine  FANDM 

This  subroutine  is  used  in  both  the  steady  and  unsteady  flow  calculations.  It 
is  called  by  the  MAIN  program  immediately  after  the  pressure  distribution  over  the 
airfoil  panels  are  known  so  that  it  can  perform  the  simple  integration  of  pressure  in  the 
appropriare  directions  to  give  the  aerodynamic  force  and  moment  coefficients  of  lift, 
drag  and  pitching  moment  about  the  leading  edge  according  to  Equations  2.22  through 
2.26. 

6.  Subroutine  GAUSS 

This  subroutine  is  the  standard  linear  system  solver  that  employs  the  well 
known  Gaussian  elimination  with  partial  pivoting  ;nd  ) Derates  simultaneously  on  i 
user  specified  number  of  right-hanu-sides.  it  is  called  by  the  MAIN  program  in  both 
the  steady  and  unsteady  flow  calculations.  In  order  to  use  GAUSS,  the  coefficients  of 
the  augmented  matrix  must  be  set  up  so  that  GAUSS  will  return  the  solutions 
replacing  the  corresponding  columns  of  the  augmented  matrix  that  were  initially 
occupied  by  the  right-hand-sides.  The  coefficient  set-ups  are  done  by  subroutines 
COFISH  and  COEFF  respectively  for  the  steady  and  unsteady  flow  problems. 


7.  Subroutine  INDATA 

This  subroutine  is  called  by  the  MAIN  program  to  read  in  the  first  three  sets 
of  data  cards  and  returns  to  the  MAIN  program  if  IFLAG  *  0.  Otherwise  it  continues 
to  read  in  the  fourth  data  card  as  the  NACA  number  corresponding  to  the  type  of 
airfoil  and  calculates  the  thickness  parameters  that  wiil  be  used  by  suoroutine 

NACA45. 

8.  Subroutine  INFL 

This  subroutine  is  the  generator  for  all  the  influence  coefficients  that  need  to 
be  stored  and  used  by  many  subroutines  associated  with  the  unsteady  flow  calculations. 
It  utilises  the  known  relative  geometrical  parameters  of  the  singularities  to  carry  out 
computations  based  on  Equations  2.12  through  2.15,  3.11  and  3.12.  The  MAIN 
program  calls  this  subroutine  in  even'  iteration  cycle  of  each  time  step  so  that  the  time- 
dependent  coefficients  can  be  updated  as  and  when  necessary.  Time-independent 
coefficients  are  computed  only  once  in  the  entire  unsteady  flow  solutions.  Those 
influence  coefficients  involving  the  wake  core  vortices  are  updated  in  each  time  step 
while  those  involving  the  shed  vorticity  panei  are  calculated  as  frequently  as  the 
number  of  iterations  take  to  terminate  a  converged  solution.  It,  however,  ices  not 
compute  and  store  those  influence  coefficients  needed  for  the  determination  of 
disturbance  potential  (Equation  3.20)  simply  because  they  are  used  only  once  in  each 
time  step. 

9.  Subroutine  KUTTA 

This  subroutine  is  called,  in  the  unsteady  flow  calculations,  by  the  MAIN 
program  during  every  iteration  cycle  in  each  time  step  to  invoke  the  Kutta  condition 
for  unsteady  flow  expressed  in  Equation  3.6.  It  calculates  the  tangential  velocities  at 
the  trailing  edge  panels  using  Equation  3.16  in  terms  of  the  unknown  vorticity  strength 
that  is  manipulated  and  solved  algebraicly. 

10.  Subroutine  NACA45 

This  subroutine  is  called  by  subroutine  BODY  if  the  airfoil  selected  belongs  to 
the  families  of  NACA  4-dig:ts  airroiis  or  the  N.ACA  5-digits  .urfoiis  of  ""pe  23dXX 
who  share  common  thickness  distributions  with  the  4-digits  airfoiis  having  the  same 
thickness  to  chord  ratio.  The  thickness  and  camber  distribution  data  of  these  airfoils 
are  calculated  and  returned  to  BODY. 


11.  Subroutine  PRESS 

This  subroutine  is  called  by  the  MAIN  program  to  calculate  the  pressure 
distribution  over  the  airfoil  panels  after  the  iterative  solution  for  the  unsteady  flow 
problem  has  successfully  met  the  convergence  criterion.  It  first  computes  the 
tangential  velocities  at  all  panel  control  pomts  using  Equation  3.16,  then  performs  the 
disturbance  potential  evaluation  at  the  current  time  step  according  to  Equations  3.20 
through  3.24.  Together  with  the  disturbance  potential  data  obtained  from  the  previous 
time  step,  it  calculates  the  pressure  distribution  using  Equation  3.19. 

12.  Subroutine  SETUP 

This  subroutine  sets  up  the  panel  nodal  coordinates  for  MAIN  program  by 
reading  the  4th  and  5th  data  sets  of  the  input  file  if  IFLAG  =  1  is  set.  It  skips  the  data 
reading  if  IFLAG  =  0  and  proceeds  to  set  up  the  node  distribution  and  call  subroutine 
BODY  to  calculate  the  airfoil  coordinates.  The  node  distribution  adopts  a  cosine 
formula  in  order  to  have  closely  packed  panels  towards  the  leading  and  trailing  edges 
for  improvements  in  solution  accuracy.  Regardless  of  how  the  nodal  coordinates  are 
obtained.  SETUP  determines  the  panel  slopes  and  airfoil  perimeter  length. 

13.  Subroutine  TEW AK 

This  subroutine  is  called  by  the  MAIN  program  at  every  iteration  cycle  of 
each  time  step  of  the  unsteady  flow  calculations  to  compute  the  resultant  velocity- 
components  at  the  mid  point  of  the  shed  vorticitv  panel  using  Equation  3.17  and  3.18 
These  velocity  components  are  necessary  to  ensure  the  correct  establishment  of  the 
shed  vorticity  panel  length  and  orientation  which  governs  the  successful 
implementation  of  the  iterative  solution  scheme  for  the  unsteady  flow  problems. 

14.  Subroutine  VELDIS 

This  subroutine  returns  to  the  MAIN  program  the  velocities  and  pressure 
distributions  for  steady  flow  calculation  using  Equations  2.20  and  2.21.  It  also 
performs  the  evaluation  of  the  disturbance  potential  at  the  panel  control  points. 
Though  these  disturbance  potential  data  are  not  necessary  for  steady  flow  solution. 

they  will  be  needed  in  the  first  time  step  of  tne  unsteady  flow  pressure  calculation. 

C.  INPUT  DATA  FOR  PROGRAM  U2DIIF 

Program  U2DIIF  reads  its  input  data  from  filecode  1.  An  example  of  the  input 
data  file  is  as  shown  in  Appendix  B  for  the  case  where  the  airfoil  nodal  coordinates  are 
input  by  user.  User  could  however  let  the  program  generate  the  nodal  coordinates  if 
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the  airfoil  chosen  happens  to  belong  to  the  family  of  \  ACA  4-digits  or  5-digits  of  type 
230XX.  To  do  this,  simply  change  I  FLAG  to  zero  m  the  first  item  of  the  3ra  set  of 
data  card  and  replace  the  nodal  coordinates  data  in  the  d01  and  oIJ1  sets  of  data  cards 
by  a  single  data  card  containing  only  the  particular  airfoil's  NACA  number  using 
Fermat  ilo').  Figure  4.2  contains  an  itemised  description  of  the  sequential  input 
variables. 

D.  OUTPUT  DATA  FROM  PROGRAM  U2DIIF 

Appendix  C  contains  a  sample  output  data  generated  by  using  the  input  data  set 
shown  in  Appendix  B.  Due  to  the  repetitive  nature  of  output  as  the  computation  time 
progresses,  only  data  at  selective  time  steps  are  shown.  The  output  data  file  begins  with 
writing  out  what  the  program  has  read  from  the  input  data  file  followed  by  the 
computed  nodal  coordinates  only  if  they  are  program  generated,  otherwise  proceeds  to 
write  the  computed  airfoil  perimeter  length.  The  next  set  of  output  data  are  the  steady 
flow  solution  parameters  of  distributed  source  strengths,  voracity  strength,  pressure 
and  velocity  distributions  as  well  as  the  force  and  moment  coefficients.  The  output  data 
terminates  at  this  point  unless  unsteady  flow  solution  is  required.  It  -hen  prints,  for 
each  time  step,  the  unsteady  Tow  solution  parameters  similar  to  the  previous  output 
for  steady  Tow  with  additional  information  pertaining  to  the  rigid  body  motions  and 
trailing  wake  vortices  data.  An  explanation  of  the  output  variable  names  are  listed  in 
Figure  4.3.  All  output  parameters  are  ncn-dimensional  quantities. 


Data  Set  =1  Format  (15)  -  I  data  card 

{TITLE  -  Number  of  utie  cards  to  be  used  in  Data  Set  =2. 

Data  Set  =2  Format  (20A4)  -  [TITLE  data  cards 

TITLE  -  Headings  to  oe  printed  on  output  for  case  run  identification. 

Data  Set  #3  Format  (315)  -  L  data  card 

I  FLAG  -  0  if  airfoil  is  NACA  XXXX  or  230XX. 

-  1  otherwise. 

NLOWER  -  Number  of  panels  used  on  airfoil  lower  surface. 

NUPPER  -  Number  of  paneis  used  on  airfoil  upper  surface  (need  not 
be  the  same  as  NLOWER). 

Data  Set  =4  Format  (6F10.6)  if  IFLAG=  l  -  variable  data  cards 

X(I)  -  x-nodal  coordinates  (divided  by  the  chord  length,  c).  A  total 

of  n  +  1  nodal  points  divided  into  6  points  per  data  card. 

Data  Set  =5  Format  (6F10.6)  -  variable  data  cards. 

VII )  -  v-nodai  coordinates  -.divided  bv  ci  corresponding  to  the 

Data  Set  =4  if  IF  LAG  =  1. 

Data  Set  =6  Format('TF10.6)  -  i  data  card 

ALPI  -  Initial  angle  of  attack  (AOA)  in  deg. 

DAL?  -  Increment  in  AOA  ;n  deg  for  non-oscillatorv  motions. 

-  Maximum  ammitude  oi  AOA  change  in  deg  lor  rotational 
harmonic  motions. 

TCON  -  Non-dimensional  rise  time  IV  X;t/c)  of  AOA  for  motion 
involving  modmed-ramp  change  in  AOA. 

FREQ  -  Non-dimensional  oscillation  frequency  (coc/V^)  for 
harmonic  motions. 

PIVOT  -  The  length  from  the  pivot  point  to  the  leading  edge  divided 
by  c  (postive  aft)  tor ‘rotational  motions. 

L'GUST  -  Magnitude  of  gust  velocity  (divided  by  V^)  along  Vqq. 

VGUST  -  Magnitude  of  gust  velocity  (divided  by  V^)  perpendicular  to  Vqq 

Data  Set  #7  Format  (3F10.3)  -  1  data  card 

DELHX  -  Amplitude  of  chordwise  translational  oscillation  divided  bv  c. 
(positive  forward). 

DELHY  -  Amplitude  of  transverse  translational  oscillation  divided  by  c. 
(positive  downward). 

PHASE  -  Phase  angle  in  deg  between  the  chordwise  and  transverse 
translational  oscillation  with  the  latter  as  reference. 

Data  Set  £8  Format  (4F10.3)  -  1  data  card 

TF  -  Final  non-dimensional  time  to  terminate  unsteady  flow  solution. 

QTS  -  Starting  time-steo  size  for  non-oscilatcry  motions  :fTADJ  =  ').0. 

-  Number  of  computation  steps  per  cvcie  for  harmonic  motions. 

-  Baseline  time-step  size  for  all  motions  if  TADJ*0.0. 

TOL  -  Tolerance  criterion  for  checking  the  convergence  between 

successive  iterations  of  (Lw)k  and  (Vw)k 

TADJ  -  Factor  by  which  DTS  will  be  adjusted. 


Figure  4.2  List  of  Input  Variables. 


TK 

TKM1 

ALPHA(T) 

OMEGA(T) 

U(T) 

V(T) 

NITR 

VXW 

VYW 

Wake 

THETA 

GAMMA 

J 

X(J) 

Q(J) 

CP(J) 

V(J) 

CD 

CL 

CM 

M 

X(M) 

Y(M) 

CIRC 


-  Time  step  tk. 

-  Time  step  tk-1. 

-  Angle  of  attack  at  time  tk. 

-  Rotational  velocity  (positive  counter  clockwise)  at  time  L. 

-  Chordwise  translation  velocity  (postive  forward)  at  time  tk. 

-  Transverse  translational  velocity  (positive  downward)  at  time  tk. 

-  Iteration  number. 

-  Iterative  solution  of(U.v)k. 

-  Iterative  solution  of  (VI, 

-  Iterative  solution  of  shed  vorticity  panel  length  A. 

-  Iterative  solution  of  shed  vorticity  panel  orientation  0,... 

-  Iterative  solution  of  the  strength  of  vorticity  distribution. 

*  Panel  number. 

-  x-coordinate  of  the  mid  point  of  j*  panel. 

-  Strength  of  source  distribution  on  the  j1*1  panel. 

-  Pressure  coefficient  at  the  mid  point  of  j*  panel. 

-  Total  tangential  velocity  at  the  mid  point  of  j*  panel 
referenced  to  the  airfoil-fixed  coordinate  system. 

-  Drag  coefficient. 

-  Lift  coefficient. 

-  Pitching  moment  coefficient  about  the  leading  edge. 

-  Trailing  wake  core  vortex  number. 

-  x-coordinate  of  the  center  of  m^  core  vortex. 

-  y-coordinate  of  the  center  of  m^  core  vortex. 

-  Circulation  strength  of  the  m*  core  vortex. 


Figure  4.3  List  of  Output  Variables. 


V.  RESULTS  AND  DISCUSSIONS  ON  CASE-RUNS 


This  Chapter  presents  the  results  of  numerous  case-runs  of  U2DIIF  code  for  the 
purpose  of  verifying  the  code.  The  various  airfoils  used  in  the  case-runs  are  deliberately 
chosen  to  be  the  same  as  those  airfoils  where  direct  comparison  of  results  can  be  made 
with  either  theoretical  analyses  and/or  experimental  data  available  in  the  literature. 

A.  STEP  CHANGE  IN  ANGLE-OF-ATTACX 

1 .  Case- Run  Definitions 

Consider  an  airfoil  initially  at  zero  angle  of  attack  to  the  free  stream  that 
undergoes  a  step  change  in  angle  of  attack  at  t~.  The  resulting  flow  should  then 
provide  the  time-dependent  information  on  the  build-up  of  aerodynamic  forces  and 
moments  on  the  airfoil  resembling  the  classical  results  of  Wagner  [Ref.  5]  calculated 
based  on  linearised  theory.  Although  Wagner  prescribed  a  slightly  different  initial 
condition  in  that  the  airfoil  is  initially  at  rest  and  impulsively  started  at  an  angle  of 
attack  and  velocity  V^,  the  difference  is  insignificant,  especially  for  a  symmetrical 
airfoil.  This  is  because  the  seemingly  different  initial  conditions  when  translated  into 
the  mathematical  model  means  that  the  step  change  in  AOA  uses  non  zero  initial 
circulation  Tg  at  tQ  with  non  zero  initial  disturbance  potential  at  infinity  if  the  airfoil  is 
cambered.  For  a  symmetric  airfoil,  these  initial  values  are  all  zeroes  and  therefore 
mathematically  would  be  the  same  as  the  initial  conditions  prescribed  by  Wagner. 

2.  Results  and  Discussions 

a.  Von  Mises  8.4%  Thick  Symmetrical  Airfoil 

A  8.4%  thick  symmetrical  Von  Mises  airfoil  is  used  for  this  case-run  where 
the  airfoil  performs  a  0.1  rad  (or  5.73°)  step  change  in  AOA.  Figure  5.1  illustrates  the 
changes  in  the  pressure  distributions  over  the  airfoil  at  time  instances  corresponding  to 
the  airfoil  having  traveled  distances,  in  terms  of  chord  length,  of  0.2,  0.5.  1.0.  2.0  and 
tc.  The  associated  trailing  wake  patterns  at  these  time  instances  i  less  :=  -ire  shown 
in  Figure  5.2.  The  time-dependent  build-up  of  aerodynamic  coefficients  of  lift,  drag, 
pitching  moment  and  the  circulation  strength  over  a  computation  period  of  two 
traveled  chord  length  are  shown  in  Figure  5.3.  Notice  that  the  lift,  pitching  moment 
and  circulation  results  are  normalised  by  the  respective  steady  state  values  at  the  same 
AOA.  The  apparently  large  initial  loading  on  the  airfoil  shown  in  Figure  5.3  correlates 


consistently  with  the  results  of  the  Piston  Theory  of  [Ref.  61  which  predicts  the  starting 
load  on  an  arbitrary  wing  to  be 

CL  =  4a  M 

‘'starting 

where  M  is  the  Mach  number.  In  the  case  of  an  incompressible  How  (M  =  0),  the  initial 
loading  would  be  infinitely  large.  The  same  large  initial  loading  was  obtained  by  Kim 
and  Mook  [Ref.  7]  who  used  continuous  vorticities  as  panel  singularities  instead  of  our 
source  and  voracity  approach.  Perhaps  what  remains  most  surprising  is  that  the  work, 
of  Basu  and  Hancock.  LRef.  3]  did  not  predict  this  initial  loading,  although  they  used 
the  same  singularity  distributions  as  U2DIIF  code.  The  initial  large  loading  in  lift  and 
pitching  moment  decreases  rapidly  over  a  short  time  span,  whereby  the  airfoil  traveled 
approximately  one-tenth  chord  length,  before  rising  in  a  manner  parallel  to  the  Wagner 
Function.  The  drag,  however,  continues  to  decrease  monotonically  after  the  initial 
sharp  fall.  The  circulation  rises,  as  continuous  shedding  of  vorticity  takes  place,  slowly 
from  the  initial  condition  of  zero  to  the  asymptotic  steady  state  value  as  time 
approaches  ».  These  results,  disregarding  the  initial  large  loading  associated  with 
incompressible  flow,  are  in  close  agreement  with  the  results  of  [Refs.  3,4,7], 
b.  Thickness  Effects  on  the  Wagner  Function 

In  order  to  more  closely  correlate  the  results  of  U2DIIF  code  to  the 
theoretical  prediction  of  Wagner,  we  performed  the  step  AOA  change  calculations  for  a 
very  thin  (1%  thickness)  NACA  4-digit  symmetrical  airfoil  which  in  reality  should 
represent  a  flat  plate.  The  results  are  plotted  as  shown  in  Figure  5.4.  Shown  also  on  the 
same  Figure  are  the  results  of  the  8.4%  thick  Von  Mises  airfoil  and  a  25.5%  thick 
symmetric  Joukowski  airfoil.  The  initial  loading  falls  off  less  rapidly  for  the  case  of  the 
simulated  flat  plate  as  compared  to  other  thick  airfoils  but  the  subsequent  rise  in  lift 
follows  very  closely  the  Wagner  Function. 
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Figure  5.2  Trailing  Wake  Patterns  at  Various  Time  Instances 
Resulting  from  a  0.1  rad  Step  Change  in  AOA  for  a 
8.4%  Thick  Symmetric  Von  Mises  Airfoil. 
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Wagner  Function 
1%  Thick  NACA-OOOl  Airfoil 
S.4°o  Thick  Von  Mises  Airfoil 
25.5%  Thick  Joukowski  Airfoil 


Figure  5.4  Time-Dependent  Lift  Resulting  from 
Step  Change  in  AOA  for  Airfoils 
of  Various  Thicknesses. 
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B.  MODIFIED  RAMP  CHANGE  IN  ANGLE-OF-ATTACK 
1.  Case-Run  Definitions 

The  case  of  a  step  change  in  AOA  can  be  considered  as  an  useful  check  for 
U2DIIF  code  since  a  handful  of  results  from  other  theoretical  analyses  are  available. 
However,  a  step  change  ;n  AOA  is  practically  not  realisable  since  all  motions,  short  of 
having  infinite  velocities,  take  piace  over  a  Unite  time  span.  The  step  change  in  AOA 
that  is  practically  possible  is  in  fact  some  form  of  ramp  rise  over  a  short  time  span  with 
large  velocity.  Even  so.  due  to  the  inertia  of  the  airfoil,  an  exact  ramp  rise  in  AOA 
does  not  physically  describe  the  actual  motion  of  the  airfoil  since  finite  time  is  also 
involved  before  the  airfoil  couid  buiid  up  its  ramp  velocity.  Same  argument  holds  at  the 
end  of  the  ramp  rise  before  the  airfoil  could  stop  at  the  final  value  of  AOA.  Therefore 
a  so  called  modified  ramp,  with  some  form  of  rounding  at  the  two  ends  of  a  ramp,  is 
mere  likely  to  describe  anything  close  to  what  is  physically  achievable.  The  theoretical 
work  of  Homentcovschi  in  [Ref.  S]  considers  the  case  of  a  flat  plate  that  moves  with 
constant  velocity  and  changes  the  incidence  about  the  mid  chord,  through  a  particular 
•amp  fashion,  described  mathematically  as. 

'0  for  t  <  0, 

a(t)  =<  6a  (3  -  2t;t)  t2;t2  for  0  ^  t  <t 
,6a  for  t  >  t 


where  6a  is  the  magnitude  of  the  AOA  change  and  T  is  the  rise  time  for  the  AOA  to 
reach  its  final  value.  This  particular  function,  plotted  as  shown  in  Figure  5.5,  does  in 
fact  describe  such  a  modified  ramp. 

2.  Results  and  Discussions 
a.  Flat- Plate  Case- Run 

Since  the  results  of  [Ref.  8]  serves  as  another  excellent  source  for  the 
verification  of  U2DIIF  code,  the  obvious  thing  to  do  is  to  use  U2DIIF  to  compute  for 

'he  case  of  i  ilat  plate,  again  simulated  by  die  !°*  thick  NACA-0001  airfoil,  ixecutmg 
'ms  modified  ramp  nse  of  A I  rad  AOA  over  a  r.se  time  of  1.5  chord  length.  This  rise 

time  is  chosen  simply  to  facilitate  a  direct  comparison  of  results  to  [Ref.  8]  which  used 
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a  rise  time  of  3  half-chord  lengths.  The  results  of  computation  are  shown  in  Figure  5.6 
and  5.7.  Figure  5.6(a)  takes  a  close  look  at  the  build  up  of  lift  during  the  initial  period 
when  the  airfoil  moves  a  distance  of  six  -herd  lengths.  The  lift  initially  rises  to  about 
82%  and  then  decreases  to  about  66 %  of  the  steady  state  value  during  the  transient 
rise  time.  Thereafter  it  increases  monotonically  in  a  manner  parallel  to  the  Wagner 
Function.  Figure  5.6(b)  is  a  zoom  view  of  the  rather  slow  convergence  of  lift  to  the 
steady  state  value.  It  takes  the  airfoil  to  cover  a  distance  of  around  50  chord  length 
before  the  lift  builds  up  to  almost  99%  of  the  steady  state  value.  The  same  results  were 
obtained  in  the  theoretical  analysis  of  [Ref.  8J.  Figure  5.7  shows  a  collection  of  the 
time-dependent  aerodynamic  parameters  resulting  from  this  particular  case-run. 
b .  Thickness  Effects 

The  same  modified  ramp  function  is  used  on  the  8.4%  thick  Von  Mises 
airfoil.  The  resulting  lift-history  plotted  in  Figure  5.8  shows  a  lower  peak  value  of  lift 
during  the  transient  AOA  rise  as  compared  to  the  case  of  a  tlate  piate  though  a  similar 
trend  of  lift  rise  is  obtained.  Figures  5.9  and  5.10  show  the  results  of  pressure 
distributions  and  trailing  wake  patterns  at  various  time  instances.  One  could  directly 
compare  these  Figures  to  the  corresponding  Figures  arising  from  the  step  AOA  change 
calculations  and  see  the  remarkable  differences  in  transient  characteristics  as  a  result  of 
varying  the  prescribed  motions.  Incidentally,  one  should  realise  that  the  non- 
dimensional  rise  time  of  1.5  chord  length  is  a  deceivingly  large  number.  In  fact,  when 
one  converts  this  to  the  real  time  for  an  airfoil  of  10  ft  chord  length  moving  at  a  low 
Mach  number  of  0.2,  the  rise  time  is  indeed  only  of  the  order  of  0.06  sec.  which  for 
practical  purpose  is  close  enough  to  a  step.  Neverthless,  the  transient  part  of  the  lift 
response  is  entirely  governed  by  how  one  prescribes  the  transient  motion. 


Figure  5,6  Normalised  Lift  C^Cg  Resulting  from  a 
Modified  Ramp  AOA  (6a  =  0.1  raS,  t  =  1.5)  about 
the  Mid  Chord  of  a  NACA-0001  Airfoil. 


Figure  5.3  Normalised  Lift  Cg  Resulting  from  a 
Modified  Ramp  AOA  (6a  =  0.1  rad.  t  =  1.5)  about 
the  Mid  Chord  of  a  Miscs  S.4%  Thick  airfoil. 


C.  TRANSLATIONAL  HARMONIC  OSCILLATION 


I.  Case- Run  Definitions 

Although  the  U2DIIF  code  is  capabie  of  computing  unsteady  flow  <c;ur:on 
for  any  general  translational  harmonic  motion  described  by  a  chordwise  and  a 
transverse  components  bearing  a  given  phase  relationship. 


hy(t)  =  5hy  sin(o)t) 
hx(ti  =  6hx  sin(<ot  +  X) 


where  a)  is  the  oscillation  frequency,  X  is  the  phase  angle  between  the  two  oscillation 

components  and  Sh  &  Sh  are  the  magnitudes  of  chordwise  and  transverse  oscillations 
x  y 

respectively.  The  case-run  to  be  presented  in  this  section  selects  the  motion  to  consist 
of  only  the  transverse  component,  i.e.  the  heaving  or  plunging  motion.  A  \AC\-d0l5 
airfoil  is  chosen  for  the  case-run.  The  airfoil  is  initially  at  zero  AOA  with  the 
freestream  Vqq  and  performs  the  plunging  oscillation  at  an  amplitude  of  6hv  and  a 
reduced  frequency  of  toe:  V x 
2.  Results  and  Discussions 

Figures  5.11  and  5.12  present  the  results  of  an  airfoil  executing  a  plunging 
motion  at  an  amplitude  of  0.018c  but  with  two  different  reduced  frequencies  of  4.3  and 
17.0  respectively.  These  numbers  are  chosen  to  coincide  with  those  numbers  used  in 
[Ref.  4],  Excellent  correlations  are  obtained.  Notice  from  these  Figures  that  the 
oscillation  frequency  has  a  great  influence  on  the  magnitudes  of  the  aerodynamic 
parameters  due  to  the  formation  of  significantly  different  trailing  wake  patterns  for  the 
same  oscillation  amplitude.  Also  to  note  is  that  the  width  of  the  resulting  trailing  wake 
is  much  larger  than  the  amplitude  of  the  oscillation,  reinforcing  the  fact  that  the 
unsteady  flow  is  strongly  governed  by  the  shed  vorticity  in  the  trailing  wake.  The  lift 
and  pitching  moment  oscillate  at  the  same  frequency  as  the  airfoil  motion  but  slightly 
out  of  phase,  the  phase  differences  vary  with  the  oscillation  frequency.  The  drag  is 


however  osciiiatina  .it  about  twice  the  frcauencv  tf  'he  airfoil  motion  with  t  negative 


mean  value,  indicating  that  the  plunging  action  indeed  generates  seme  propulsive 

thrust.  The  same  conclusion  was  arrived  at  in  the  experimental  work  of  Halfman 
[Ref.  9]  using  a  symmetrical  NACA-0012  airfoil. 


(a)  Position  of  Airfoil  (Positive  Downward) 


■10  2.3  J.o  +o 

.cot;  (2;r; 

(b)  Time  History  of  Cj,  Cm  and  Cd  over  2  Cycles. 


ure  5.11  Harmonic  Plunging  Motion  of  a  NACA-0015  Airfoil 
0hv  =  0.018c.  Wc  Vx  =  4.3. 
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Figure  5.11  .  (cont'd.) 
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(a)  Position  of  Airfoil  (Positive  Downward) 


u)u  t  :*> 

(b)  Time  History  of  C^,  Cm  and  Cd  over  2  Cycles. 


Figure  5.12  Harmonic  Plunging  Motion  of  a  NACA-0015  Airfoil 

oh  =  0.0 1  Sc,  cdc A'-*-  =  17.0. 

> 


D.  ROTATIONAL  HARMONIC  OSCILLATION 

1.  Case- Run  Definitions 

The  case  of  an  airfoil  oscillating  harmonically  in  pitch  will  be  uniquely  defined 
only  if  a  pivot  point  is  prescribed.  If  a  fixed  pivot  point  is  used  in  the  calculations,  take 
for  example  the  leading  edge,  any  pitching  motion  about  a  pivot  point  other  than  the 
leading  edge  would  need  to  be  described  as  composed  of  a  pitching  and  a  translation 
motions  about  the  leading  edge.  Program  U2DIIF  handles  such  conversion 
automatically  without  having  the  user  to  figure  out  the  combined  motion.  This  applies 
also  to  the  case  of  modified  ramp  rise  in  AOA.  The  harmonic  pitching  oscillation  is 
described  by, 

a(t)  =  8a  sin(tot) 

where  8a  and  a>  are  the  amplitude  and  frequency  of  oscillation  respectively. 

2.  Results  and  Discussions 

The  results  for  the  case  of  the  S.4%  thick  Von  Mises  symmetric  airfoil 
oscillating  at  an  amplitude  of  0.01  rad  (or  0.573°)  at  a  rather  high  reduced  frequency 
of  toe.  =  20.0  about  the  leading  edge,  are  shown  in  Figure  5.13.  The  lift,  drag  and 
pitching  moment  time-history  as  well  as  the  trailing  wake  patterns  are  very  much 
similar  to  the  case  of  a  plunging  airfoil  at  frequency  of  the  same  order  of  magnitude. 
The  differences  are  clearly  in  the  magnitudes  and  phase  angles.  These  results  check 
closeiy  to  those  of  [Ref.  3].  Figure  5.14  shows  the  results  of  the  same  Mises  airfoil 
performing  another  harmonic  oscillation  at  a  lower  reduced  frequency  of  0.8  and  a 
large  amplitude  of  0.3973  rad  about  a  pivot  point  0.5  chord  length  ahead  of  the  leading 
edge.  [Ref.  4]  conducted  the  same  analysis  for  this  pitch  oscillation  although  the 
reason  for  using  such  a  high  amplitude  of  almost  23°  was  not  clear.  It  is  envisaged 
that  such  high  amplitude  may  result  in  flow  separation.  Nevertheless,  the  case-run  is 
carried  out  assuming  validity  of  attached  flow  for  the  sake  of  comparing  the  results. 
P-irhaos  m  inherent  advantage,  in  the  light  of  V2DIIF  code  verification,  with  the  use 
:f  .ugh  amplitude  m  tins  case-run  is  that  a  discrepancy.  .1  my,  would  mow  :p 
significantly.  Due  to  the  use  of  different  computation  time-step  size  the  pressure 
distributions  on  the  airfoil,  shown  in  Figure  5.14,  do  not  correspond  one-to-one  at 
exactly  the  same  angular  positions  as  those  presented  in  [Ref.  4],  However,  the  angular 
positions  are  matched  to  within  0.001°,  0.05°  and  0.8°  respectively  for  the  three 
pressure  distributions  shown.  They  correlate  very  well  to  the  results  of  [Ref.  4). 


Figure  5.13  Harmonic  Pitching  Motion  about  the  Leading  Edge 
of  a  8.4%  Thick  Von  Mises  Airfoil 
6a  =  0.01  rad,  wc.Vqq  =  20.0. 


Figure  5.14  Harmonic  Pitching  Motion  about  a  Pivot  0.5c  ahead  of 
the  Leading  Edge  of  a  8.4%  Thick  Von  Mises  Airfoil 
5a  =  0.3973  rad,  coc/V^  =  0.8. 


e  Ohoravr.se  Station  as  Fraction  :t'  Chord) 

(c)  Pressure  Distribution  at  a  =  0.1876  rad. 


E.  SHARP  EDGE  GUST  FIELD  PENETRATION 


1.  Case- Run  Definitions 

The  ease  of  an  airfoii  penetrating  a  sharp  edge  gust  field  can  be  computed 
using  the  U2DIIF  code  by  specifying  the  components  of  gust  velocity  along  and 
perpendicular  to  the  freestream  Vx  and  the  angle  of  attack  of  the  airfoil.  The  resuits 
that  are  presented  here  consider  the  case  of  airfoil  penetrating  a  sharp  edge  vertical 
gust  at  zero  AOA,  in  view  of  generating  information  on  the  time-dependent  lift 
resembling  the  classical  results  of  Klissner  [Ref.  10]  based  on  linearised  theory.  The 
gust  front  is  taken  to  be  at  the  leading  edge  at  tQ  with  only  the  transverse  (vertical) 
component  of  0.25V  -q . 

2.  Results  and  Discussions 

Figures  5.15  and  5.16  demonstrate  the  variation  of  pressure  distributions  and 
trading  wake  patterns  respectively  during  and  shortly  after  the  gust  front  moves  past 
the  airfoil.  It  is  interesting  to  see  that  the  resulting  wake  pattern  after  the  entire  airfoil 
is  submerged  in  the  gust  field  is  as  if  being  split  by  the  gust  front  into  two  portions 
curiing  in  opposite  directions.  Ref.  31  predicted  a  similar  behaviour  for  the  case  of  an 
nndeformed  gust  front.  Due  to  the  use  of  the  modified  flow  tangency  condition  to 
handle  the  situation  when  the  gust  front  happens  to  fall  in  between  two  nodal  points, 
the  pressure  distributions,  predicted  by  U2DIIF  code,  lie  in  between  the  results  of  the 
undeformed  and  deformed  gust  front  models  used  in  [Ref.  3].  We  therefore  conclude 
that  this  modified  flow  tangency  condition  produces  sufficiently  accurate  results 
without  adding  complication  in  deformed  gust  field  modeling  which  in  turns  limits  the 
application  to  only  sharp  edge  gusts.  The  present  method  would  therefore  preserve  the 
generality  for  extending  calculations  to  other  types  of  gust  front.  Another  comparison 
is  made,  as  showm  in  Figure  5.17,  by  plotting  the  build-up  of  lift  as  a  function  of 
distance  traveled  by  the  airfoil  in  chord  lengths.  Shown  in  the  same  Figure  are  the 
Klissner  Function  and  the  results  obtained  from  another  case-run  using  a  25.5° o  thick 
Joukowski  airfoil. 


x  c  (Chordwise  Station  as  Fraction  of  Chord)  • 


(a)  tV'gQ,  c  =  0.25 


Figure  5.]5  Pressure  Distributions  at  Various  Time  Instances 
from  a  8.4%  Thick  Von  Mises  Airfoil  Penetrating  a 
Vertical  Sharp  Edge  Gust  of  0.25%^. 
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(c)  tVgo/C  =  1.0 


Figure  5.15  .  (cont'd.) 


Figure  5.16  Trailing  Wake  Patterns  at  Various  Time  Instances  Resulting 
from  a  S.4%  Thick  Von  Mises  Airfoil  Penetrating  a 
Vertical  Sharp  Edge  Gust  of  0.25V X). 
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(d)  tVgo/c  *  2.0 
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Figure  5.16  .(corn'd.) 
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Kiissner  Function 

8.4%  Thick  Von  Mises  Airfoil 

25.5%  Thick  Joukowski  Airfoil 


Jl 


Figure  5.17  Normalised  Time- Dependent  Lift  Cf  Cf  due  to 
Airfoils  ot  Various  Thicknesses  Penetrating  00 
a  Vertical  Sharp  Edge  Gust  of  0.25VX. 
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VI.  CONCLUDING  REMARKS 


A.  GENERAL  COMMENTS 

The  U2DIIF  computer  code  has  been  developed  for  the  purpose  of 
demonstrating  the  successful  extension  of  the  well  known  panel  methods,  which  have 
been  used  extensively  for  steady  flow  problems,  into  a  powerful  numerical  tool  for 
solving  the  unsteady  flow  problems.  The  mathematical  modeling  of  the  various  types 
of  unsteady  flows  has  been  done  with  the  goal  of  preserving  the  generality  of  the 
methods.  The  intention  is  to  present  a  method  that  has  the  minimum  inherent 
limitations  and  restrictions  so  that  its  usage  for  future  applications  and  developments 
could  remain  appealing. 

The  validation  of  the  U2DIIF  code  has  been  done  through  the  various  case-runs 
of  the  numerous  types  of  unsteady  flow  problems.  The  results  of  each  case-run  has 
been  shown  to  be  well  correlated  to  the  results  obtainable  from  the  literature  .n  'he 
form  of  theoretical  analyses,  numerical  calculations  based  on  different  variants  of  panel 
singularities  and  in  cases  where  limited  experimental  data  are  available.  The  ability  of 
those  case-runs  using  an  airfoil  as  thin  as  1%  thickness  to  produce  results  that 
correlate  accurately  to  the  theoretical  flat  plate  results  is  perhaps  a  remarkable 
robustness  possessed  by  the  present  unsteady  flow  solution  methods. 

B.  ENHANCING  U2DIIF  PROGRAM'S  CAPABILITY 

It  has  been  noted  in  Chapter  IV  that  the  current  U2DIIF  code  limits  the  total 
number  of  panels  to  200  and  the  total  computation  time  steps  to  200  also.  These  limits 
are  not  at  all  rigid  and  can  be  easily  increased  if  the  computer  storage  space  is  not 
critical.  A  point  to  note  is  that  the  computation  time  will  grow  rapidly  as  these  limits 
are  raised.  The  current  linear  system  solver,  the  Gaussian  elimination  algorithm,  used 
in  the  code  must  be  concurrently  improved  upon  to  efficiently  reduce  the  computation 
■.me  required.  for  me  iterations  .a  each  time  >tep.  A  close  examination  of  :ne  natrix 
Equation  3.15,  where  the  linear  system  solver  is  needed  in  every  iteration,  reveals  that 
the  coefficients  of  the  left-hand-side  matrix  [  A  ]  are  time-independent  constants. 
Therefore  the  Gaussian  elimination  algorithm  need  only  be  done  once  for  the  entire 
unsteady  flow  calculations  as  far  as  the  left-hand-side  matrix  is  concerned.  One  could 
then  perform,  for  each  iteration,  the  manipulation  of  the  two  right-hand-sides 


according  to  the  steps  taken  for  the  reduction  of  the  left-hand-side.  This  should  cut 
down  the  computation  time  significantly  against  the  current  method  of  manipulating 
both  the  left-  and  right-hand-sides  simultaneously  in  each  iteration.  The  savings  in 
computation  time  was  not  pursued  in  the  development  of  U2DIIF  code  because  the 
prime  concern  was  to  demonstrate  that  the  basic  iterative  solution  scheme  works  for 
the  unsteady  flow  problems. 

Another  improvement  that  one  may  consider  is  to  enable  the  code  to  be 
continuable  from  a  time  step  where  previous  calculations  were  terminated.  One  sees 
this  requirement  necessary  not  only  in  the  case  of  premature  termination  of 
computation  due  to  some  unforeseen  circumstances,  but  also  if  one  needs  to  proiong 
the  computation  time.  Certainly  with  the  current  code  structure,  one  has  no  choice  but 
to  perform  the  calculation  right  from  the  beginning. 

A  farther  extension  of  U2DIIF  code  to  the  computation  of  unsteady  flow 
problems  involving  multiple  airfoils  may  be  tvorth  pursuing.  Other  research  works  that 
could  be  done  based  on  U2DIIF  code  are  in  the  area  of  incorporating  more  variety  of 
rigid  body  motions  into  the  code  either  in  the  form  of  closed-form  equations  or 
tabulated  time  history  of  the  positions  and  rates  of  motions.  It  is  important  that  one 
should  use  as  close  as  possible,  in  the  code,  a  rigid  body  motion  that  describes  the 
physical  motion  before  generating  any  numerical  unsteady  flow  results  for  meaningful 
comparison  to  test  data.  This  fact  has  been  well  illustrated  and  emphasized  in  the 
comparison  of  results  of  case-runs  involving  a  step  change  to  that  of  a  ramp  change  in 
AOA  with  a  fast  rate  which  one  could  regard  as  a  step  in  reality.  However,  remarkable 
difference  in  transient  aerodynamics  has  been  shown. 


APPENDIX  A 

U2DIIF  PROGRAM  LISTINGS 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  PROGRAM  U2DIIF 

c 

C  UNSTEADY  MOTION  OF  A  TWO-DIMENSIONAL  AIRFOIL 

C  IN  INCOMPRESSIBLE  INVISCID  FLOW 

C  USING  PANEL  METHODS  BASED  ON  THE  HESS  &  SMITH 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON  /BOD/  IFLAG .NLOWER , MUPPER ,N0DT0T,X(202) , Y(202) , 

+  COSTHE ( 201 ) , 3INTKE (201) , SS ,NP1 ,NP2 

COMMON  / WAK/  VYW , VXW , WAKE , DT 
COMMON  /WAK2/  VYWK.VXWK 
COMMON  /SING/  Q (200) .GAMMA, QK(200) ,GAMK 

COMMON  /COM/  CV(200 )  ,XC(200)  , YC(200)  ,M,TD,CWX(200)  ,CWY(200) 
COMMON  / POT/  PHI (200 ) , ?HIX(200 ) 

COMMON  /GUST/  UG 1 2C0 )  ,  7G  ( 200 )  ,  XGF ,  UGUST ,  VG'JST 
COMMON  /EXTV/  UE(200) 

DIMENSION  XXC(200) , YYC(20Q) 

n 

C  INPUT  FROM  FILE  CODE  5  AND  SET  UP  PANEL  NODES  AND  SLOPES 

PI  =  3 . 1415925535 

WRITE  ( 6. 1003) 

1003  FORMAT  (/////,•  DATA  READ  FROM  FILE  CODE  1',//) 

CALL  INDATA 
CALL  SETUP 

READ  (1,501)  ALP I, DALP,TCON, FREQ, PIVOT, UGUST, VGUST 
WRITE  (6.501)  ALP I, D ALP, TCON, FREQ, PIVOT, UGUST, VGUST 
501  FORMAT  (7F1Q.6) 

READ  (1,501)  DELHX.DELHY, PHASE 
WRITE  (6.501)  DELHX.DELHY, PHASE 
READ  (1,501)  TF , DTS , TOL , TAD J 
WRITE  (6,501)  TF, DTS, TOL, TAD J 
IF  (IFLAG  .EQ.  0)  WRITE  (6,1005) 

1005  FORMAT  (///,'  COORDINATES  OF  AIRFOIL  NODES', 

+  /  /  ,3X,  '  X/C  ,6X  '  Y/C  ,/) 

IF  (IFLAG  .EQ.  0)  WRITE  (6,1010)  (X(I ) , Y(I ) , 1=1 ,NODTOT+l ) 

1010  FORMAT  (F10 . 6 ,F10 .6) 

WRITE  (6,1020)  SS 

1020  FORMAT (// , 1  AIRFOIL  PERIMETER  LENGTH  =  ' ,F10.6,/) 

C 

C  STEADY  FLOW  CALCULATION  AT  ALP  I 
C 

ALPHA  =  ALPI 
WRITE  (6  1030)  ALPHA 

1030  FORMAT  (//,'  STEADY  FLOW  SOLUTION  AT  ALPHA  =  ' ,F10.6,/) 

IF  < ALPHA  .GT.  90.’)  GO  TO  200 
COSALF  =  COS ‘ALPHA*?!/ 130. ) 

3 INALF  =  3 IN ( ALPHAS? I / 130 . ; 

CALL  CO F I SK < S INALF , CCS ALF ) 

CALL  GAUSS ( 1 ) 

CALL  VELDIS(SINALF, COSALF) 

CALL  FANDM(SINALF, COSALF) 

C 

C  INITIALISATION  FOR  UNSTEADY  FLOW  CALCULATION  TO  BEGIN 


HX 

HY 

HXO 

HYO 


=  0.0 
=  0.0 
=  0.0 
=  0.0 
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ooooooo 


I 

* 


1 

i 


■Si 

•«* 

i- 


•c 

■?: 


1 


a 


DHX  =0.0 

DHY  =  0.0 

TJX  =  0 . 0 

UY  =0.0 

ALP  =  ALPI 
DA  =0.0  ■ 

COSDA  =  1.0 

SIND A  =0.0 

OMEGA  =  0.0 

HGF  =0.0 

ANGLE  =  AL?I*?I/130.  +  ATAN ( VGUST/ ( 1 . +UGUST ) ) 

COSANG  =  COS (ANGLE) 

SINANG  =  SIN (ANGLE) 

DO  100  IG  =  1 ,NODTOT 
UG(IG)  =  0.0 

7G('  IG)  =0.0 
?HA  =  PHASE*?!/ 130 . 

7XW  =  COSALF 

VYW  =  SINALF 

GAMK  =  GAMMA 

T  =  0.0 

M  =0.0 

TOLD  =0.0 

RIGID  30DY  MOTIONS  OF  AIRFOIL 

IF  (FREO  .ME.  0.0)  GO  TO  1 

IF  ( DAL?  .20.  0.0)  GO  TO  2 

IF  (TCON  .NE.  0.0)  GO  TO  3 

ALPHA  =  ALPI  +  DALP 
COSALF  =  COS (ALPHA*?!/ 130. ) 

SINALF  =  5 IN ( ALPHA*? I / 1 30 . j 
DT  =  DTS 

TD  =  DTS 

GO  TO  50 

IF  (OJGUST  .EQ.  0.0)  .AND.  (VGUST  .EQ.  0.0))  GO  TO  200 

DT  =  DTS 

TD  =  DTS 

GO  TO  60 

DT  =  2.0*PI/(FREQ*DTS) 

TD  =  DT 

T  =  DT 

WRITE  (6,1051) 

FORMAT  (/////,'  ***************************************  i  / 
+  1  ***  BEGIN  UNSTEADY  FLOW  SOLUTION  **** ' . / , 

+  1  +************************************** i \ 

M  =  M  +  1 

IF  (T  .GT.  TF)  GO  TO  200 

STORE  CORE  VORTEX  COORDINATES  FOR  TIME  STEP  ADJUSTMENTS 

IF  (M  .EO.  1)  GO  TO  50 
DO  51  I  =  1 ,M-1 
XXC(I)  =  XC(I) 

YYC(I)  =  YC( I ) 

I?  (FREO  .NE.  0.0)  GO  TO  11 
IF  DAL?  .EO.  0.0)  SO  TO  22 

-  *  -  — <M  .  .la .  0.0)  'jO  -0  j  j 

STEP  CHANGE  IN  AOA 

IF  (TADJ  .NE.  0.0)  GO  TO  70 
TD  =  FLOAT (M+1)*DTS 

GO  TO  70 

MODIFIED  RAMP  CHANGE  IN  AOA 

IF  (T  .GT.  TCON)  GO  TO  34 
DAL  =  DALP  *  (3.-2 .*T/TCON)*(T/TCON)**2 
ALPHA  =  ALPI  +  DAL 


34 


70 


C 

c 

c 


COSALF 
SINALF 
DA 

COSDA 
SINDA 
OMEGA 
DHX 
DHY 
UY 

MTCON 
GO  TO 
DAL  = 
ALPHA 
COSALF 
SINALF 
DA 

COSDA 

SINDA 

OMEGA 

DHX 

DHY 

UY 

IF  (TADJ 
TD 

GO  TO  70 


COS (ALPHA*PI/180 . ) 

SIN(ALPHA*?I/130. ) 

ALPHA  -  ALP 
COS(DA*PI/180. ) 

SIN(DA*PI/180.) 

-  (DALP*PI/180. )  *  (6 . *T/ (TCON*TCON) )  *  (l.-T/TCON) 
PIVOT  *  (1. -COSDA) 

*  SINDA 
OMEGA 


-  PIVOT 
PIVOT  * 
M 


0.0 

ALPI  +  DALP 
COS (ALPHA*PI/180 . ) 
SIN(ALPHA*PI/180.) 
0.0 
1.0 
0.0 
0.0 
0.0 
0.0 
0.0 

.NE.  0.0)  GO  TO  70 
■  FLOAT (M+l -MTCON )*DTS 


22 


SHARP  EDGE  GUST  (UGUST  AND/OR  VGUST) 


=  T 

IG  =  1 ,NODTOT 
=  0.0 
=  0.0 

=  X(IG)*COSALF  +  Y(IG)*SINALF 
=  X(IG+1) *C0SALF  +  Y(IG+1)*SINALF 
.LT.  NLOWER+1)  GO  TO  120 
.LE.  XG)  GO  TO  110 


XG) 


XGP1 ) 


IF  (XGF  ,GE.  XGP1 )  GO  TO  111 
FAC  =  (XGF  -  XG)  /  (XGP1  - 
UG(IG)  =  UGUST*FAC 
VG(IG)  =  VGUST*FAC 
GO  TO  110 
UG(IG)  =  UGUST 
VG(IG)  =  VGUST 
GO  TO  110 

IF  (XGF  .LE.  XGP1 )  GO  TO  110 
IF  (XGF  .GE.  XG)  GO  TO  121 
FAC  =  (XGF  -  XGP1)/ (XG  - 
UG(IG)  =  UGUST*FAC 
VG(IG)  =  VGUST*FAC 
GO  TO  110 
UG(IG)  =  UGUST 
VG(IG)  =  VGUST 
CONTINUE 

IF  (XGF  .LE.  COSALF)  MGUST  =  M 
IF  (TADJ  .NE.  0.0)  GO  TO  70 

IF  (XGF  .GT.  COSALF)  TD  =  FLOAT (M+l -MGUST ^DTS 
GO  TO  70 

TRANSLATION  HARMONIC  OSCILLATION 

IF  (DALP  .NE.  0.0)  GO  TO  12 
HX  =  DELHX  *  SIN(FREQ*T  +  PHA) 

HY  =  DELHY  *  SIN(FREQ’I'T) 

DHX  =  HX  -  HXO 

DHY  =  HY  -  HYO 

UX  =  DELHX*FREQ*COS ( FREQ*T+PHA) 

UY  =  DELHY*FREQ*COS(FREQ*T) 

GO  TO  70 

ROTATIONAL  HARMONIC  OSCILLATION 


c 

12  DAL  =  DALP*SIN  ( FREQ*T ) 

ALPHA  =  ALP I  +  DAL 
COSALF  =  COS(ALPHA*PI/180.) 

SINALF  =  SIN  (ALPHA ''PI/ 180  . ) 

DA  =  ALPHA  -  ALP 

COSDA  =  C0S(DA*PI/180. ) 

SINDA  =  SIN(DA*PI/180.) 

OMEGA  =  -  (DALP*PI/180 . )  *  FREQ  *  COS (FREQ*T) 

UY  =  PIVOT  *  OMEGA 

DHX  =  PIVOT  *  (1. -COSDA) 

DHY  =  -  PIVOT  *  SINDA 

C 

C  TRANSFORM  CORE  VORTEX  COORDINATES  W.  R.  T.  NEW  AIRFOIL  POSITION 

C 

70  IF  (M  .EQ.  1)  GO  TO  80 
DO  90  I  =  1  ,M-1 

XC(I)  =  XXC(I)  +  CWX(I)  *  DT 

YC(I)  =  YYC(l)  +  CWY(I)  *  DT 

XCO  =  XC(I) 

YCO  =  YC(l) 

XC(I)  =  XCO*COSDA  -  YCO*SINDA  +  DHX 
90  YC(I)  =  XCO*SINDA  +  YCO*COSDA  +  DHY 
80  CONTINUE 

WRITE  (6,1001)  T,DT 

1001  FORMAT  (/////,  '  TIME  STEP  TK  =  '  ,  F10 . 5 , 10X,  1 TK  -  TKM1  =  '  ,  F10 . 5  ,  /  ) 

WRITE  (6,1004)  ALPHA, OMEGA, UX,UY 
1004  FORMAT  (/,'  ALPHA (T)  =  ',F10.6,5X,'  OMEGA (T)  =  ',F10.5 ./, 

+  1  U(T)  =  ',F10.6,5X,'  V(T)  =  1 , F1C . 6 , /// , 

+  IX,'  NITR  VXW  VYW  WAKE  THETA  GAMMA',/') 

C 

C  CALCULATE  THE  TRAILING  EDGE  WAKE  ELEMENT 
C 

NITR  =  0 

10  WAKE  =  SQRT ( VYW*VYW+VXW*VXW ) *DT 
THENP1  =  AT AN 2 (VYW , VXW ) 

COSTHE(NPl)  =  COS (THENPI) 

sinthe(npi)  =  SIN (thenpi ) 

WRITE  (6,1002)  NITR, VXW, VYW, WAKE, THENPI, GAMK 

1002  FORMAT  (l5,4F10.6,E14.6) 

X(NP2)  =  X(NP1)  +  WAKE*COSTHE (NP1 ) 

y(np2)  =  y(npi)  +  WAKE*SINTHE(NP1) 

CALL  INFL  (NITR) 

CALL  COEF  (SINALF, COSALF, OMEGA, UX,UY) 

CALL  GAUSS (2) 

CALL  KUTTA  (ALPHA, SINALF, COSALF, OMEGA, UX,UY) 

CALL  TEWAK  (SINALF, COSALF) 

TOL1  =  ABS (VYW  -  VYWK) 

TOL2  =  ABS ( VXW  -  VXWK) 

IF  ( (TOL1  .LT.  TOL)  .AND.  (TOL2  .LT.  TOL))  GO  TO  20 

VYW  =  VYWK 

VXW  =  VXWK 

IF  (NITR  .GT.  200)  STOP 

NITR  =  NITR  +  1 

GO  TO  10 

20  WRITE  (6,1011)  NITR 

1011  FORMAT  CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =  .12) 

CALL  7RESS  ( SINALF ,  COSALF ,  OMEGA ,  U X ,  'JY ) 

IF  ((UGUST  .EQ.  0.0)  .AND.  . VGU5T  .EQ.  0.0))  GO  TO  300 

CALL  FANDM  (SINANG, COSANG) 

GO  TO  400 

300  CALL  FANDM  (SINALF, COSALF) 

400  CONTINUE 
C 

C  ADJUST  TIME  STEP  (TADJ  .NE.  0.0)  IF  NECESSARY 
C 

IF  (TADJ  .EQ.  0.0)  GO  TO  95 
WRITE  (5,2001) 

2001  FORMAT  (//.  1  DO  YOU  WANT  TO  ADJUST  TIME  STEP  ?  0  -  NO,  1  -  YES 1 ) 
READ  (5, IDT 


nnnnnnnnnnn  non  onn 


IF  'IDT  .SO.  0)  00  TO  95 
DT  =  TADJ  "  DT 

T  =  TOLD  -r  OT 

WRITS  (6.1006) 

1006  FORMAT  3ACK-TRACK  COMPUTATION  AND  ADJUST  TIME-STEP',//) 

GO  TO  50 

WAKE  ELEMENT  LEAVES  TRAILING  EDGE  AS  A  CORE-VORTEX 

95  C/CM)  =  S3" ' GAMMA-GAMK ! 

XC(M)  =  a(MP1 j  -  0 . 5"WAKE"C0STHE (NPI ) 

YC(M)  =  Y (NPI )  t  0 . 5 "WAKE"S INTHE (NPI) 

CWX(M)  =  VXW 
CWY(M)  »  VYW 
WRITE  (6,1052) 

1052  FORMAT  ‘.i .  TRAILING  VORTICES  DATA’  //. 

+  4X, ’M1 ,4X. 1 X ( M ) • ,6X, 'T(M) 1 (SX, 'CIRC1 ,/) 

DO  900  1=1, M 

900  WRITE  (6,1050)  I ,XC(I) ,YC(I) ,CV(I) 

1050  FORMAT ( IS , 3F10 . a  ; 

CALL  CORVOR  (SINALF , COSALF) 

RE-INITIALISE  PARAMETERS  FOR  NEXT  TIME  STEP  CALCULATION 

DO  20  I  =  l.NCDTOT 

Q ( j- )  =  0K(*; 

f'r.Zll)  =  ?HIK(I) 

30  CONTINUE 

GAMMA  =  GAMK 

AL?  =  ALPHA 

HXO  =  HX 

HYO  =  HY 

TOLD  =  T 

DT  =  TD 

T  =  T  +  TD 

GO  TO  40 

200  STOP 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  30DY(Z,SIGN,X,Y) 

RETURN  COORDINATES  OF  POINT  ON  THE  BODY  SURFACE 

Z  =  NODE-SPACING  PARAMETER 
X,Y  =  CARTESIAN  COORDINATES 
SIGN  =  +1.  FOR  UPPER  SURFACE 
-1.  FOR  LOWER  SURFACE 


CO 


c 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  BODY ( Z, SIGN, X,Y) 

COMMON  /PAR/  NACA , TAU , EPSMAX , PTMAX 
IF  (SIGN  .LT.  0.0)  Z  =  1.  -  Z 
CALL  NACA45(Z, THICK, CAMBER, BETA) 

X  =  Z  -  SIGN*THICK*SIN ( BETA) 

Y  =  CAMBER  *  SIGN*THICK*C0S(3ETA) 

RETURN 


2ND 


SUBROUTINE  COEF  (SINALF, COSALF, OMEGA, UX,UY) 

SET  COEFFICIENTS  OF  N  EOUS  ARISING  FROM  FLOW 

TANGENCY  CONDITIONS  AT  MID  POINTS  OF  PANELS 
SOLVING  THE  N- SOURCE  STRENGTHS  IN  TERMS  OF  THE 
VORTICITY  STRENGTH  (RESULTING  IN  2  RHS) 

KUTTA  CONDITION  IS  SATISFIED  SEPARATELY  TO  OBTAIN 
THE  VORTICITY  STRENGTH 

THIS  SOLUTION  METHOD  IS  DESIRED  FOR  UNSTEADY  FLOW 


CC 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 


109 


nnnnnnnnnn 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  COEF  SINALF.  COSALF, OMEGA, *JX,UY) 

COMMON  BOD/  I FLAG . NLOWER , NUPPER , NODTCT , X ( 202 ) , Y ( 202 ) , 

-  COSTHE ( 201 ) , 3INTHE (201 ) , SS ,NP1 ,NP2 

COMMON  /COF /  Ai 201,211) ,NEQS 
COMMON  /SING/  0(200) , GAMMA, QK( 200) ,GAMK 
COMMON  /WAK/  VYW , 7XW , WAKE , DT 

COMMON  / CORV/  CV(2G0) ,XC(200) ,YC(200) ,M,TD,CCVX(200) ,CCVY(200) 
COMMON  / IMF 1 /  AAN (201, 201), 3BN (20 1,201), AYNP 1(201), BYNP 1(201) 
COMMON  /IMF 2/  CCN( 20 1,200) , COT (201 , 200 ) , CYNP1 (200) , CXNP1 ( 200 ) 
COMMON  /  GUST/  rJG ( 200  )  ,  7G ( 200 )  ,  XGF ,  UGUST ,  VGUST 

NECS  =  MODTOT 

NP1  =  NODTOT  +  1 

NP2  =  NODTOT  +  2 

C 


INITIALISE  COEFFICIENTS 


90 


C 

C 

c 


DO  9C  I  =  1, NODTOT 
DO  90  J  =  1 , MP2 
A ( I , J )  =  0.0 

SET  LHS  MATRIX  A(I,J) 


DO  120  I  =  1, NODTOT 

XMID  =  0.5  *  ;  X  ( I )  +•  X(I*1) ) 

CHID  =  0.5  *  (Y(I)  ♦  Y(Ifl; ) 

3  =0.0 

DO  110  J  =  1, NODTOT 

A( I , J )  =  AAN ( I , J ) 

3  =  3  +  3BN ( I , J ) 

110  CONTINUE 


FILL  IN  THE  RIGHT  HAND  SIDE 


C 

C 

C 


100 

140 

120 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

C  SUBROUTINE  COFISH (SINALF, COSALF) 

C 

C  SET  COEFFICIENTS  OF  LINEAR  SYSTEM  -  N+l  EQUATIONS 

I  H  EOUS  -  FLOW  TANGENCY  AT  MID  POINTS  OF  PANELS 

I  .ECU  -  a ATT A  CONDITION  AT  TRAILING  EDGE  PANELS 

C  THIS  SOLUTION  METHOD  IS  EFFECTIVE  FOR  STEADY  FLOW,  MO 

C  ITERATION  IS  REQUIRED,  N-SOURCE  STRENGTHS  AND  1 

C  VORTICITY  STRENGTH  ARE  SOLVED  SIMULTANEOUSLY 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  COFISH (SINALF, COSALF) 

COMMON  /BOD/  IFLAG, NLOWER, NUPPER. NODTOT, X(202) ,Y(202) , 

+  COSTHE (201) ,SINTHE(201) ,SS,NP1 ,NP2 

COMMON  /COF/  A(201,211), KUTTA 
COMMON  /NUM/  PI,PI2INV 
KUTTA  =  NODTOT  +  1 
C 


A  :,NP1  =  - 

A(I,N?2)  =  - 
+  +  SINTHE(I) 
+  -  COSTHE (I) 


ADD  CORE  VORTEX  CONTRIBUTION 

IF  (M  .EQ.  1)  GO  TO  140 
MM1  =  M  -  1 
DO  100  N  =  1 ,MM1 

A(I ,NP2)  =  A(I ,NP2)  -  CCN(I ,N)*CV (N) 
CONTINUE 
CONTINUE 
CONTINUE 
RETURN 
END 


3  +  3BM ( I , NP1 ) *SS/WAKE 
3BN ( I , NP1 ) *GAMMA * S S / WAKE 

*  ( (1 .+UG(I) )*COSALF-VG(I)*SINALF+UX) 

*  ( ( 1 . +UG ( I ) ) *S INALF+VG ( I ) *COSALF+UY ) 

+  OMEGA* (YMID*SINTHE ( I )  +  XMID*COSTHE(I) ) 


onot  iminnnn 


c 

c 


90 


C 

C 

C 


C 

c 

c 


INITIALISE  COEFFICIENTS 

DO  90  J  =  i , KUTTA 
A ( KUTTA , J )  =  0.0 

SET  VN  =  0  AT  MID-POINT  OF  I-TH  PANEL 


DO  120  I  =  1 , NODTOT 
XMID  =  .5* 

YMID  =  .5* 

A( I, KUTTA)  =  0.0 


f  W  4  V  4 

m : 


FIND  CONTRIBUTION  OF  J-TH  PANEL 


100 


C 

c 

c 

c 


110 


DO  110  J  =  1, NODTOT 
FLOG  =0.0 
FT AN  =  PI 

IF  (J  .EQ.  I)  GO  TO  100 
DXJ  =  XMID  -  X(J) 

DXJP  =  XMID  -  X(J+1) 

DYJ  =  YMID  -  Y(J) 

DYJP  =  YMID  -  Y( J+l) 

FLOG  =  .5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
FT AN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJ?*DYJ) 

CTIMTJ  =  COSTHE ( I ) *C0STHE ( J )  +  SINTHE ( I ) *SINTHE ( J ) 

5TIMTJ  =  S INTHE ( I ) *COSTHE ( J )  -  COSTHE(I)*SINTHE(J) 

A(I , J)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 

B  =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 

A(I , KUTTA)  =  A(I, KUTTA)  +  B 

IF  ((I  .GT.  1)  .AND.  (I  .LT.  NODTOT) )GO  TO  110 

IF  I-TH  PANEL  TOUCHES  TRAILING  EDGE, 

ADD  CONTRIBUTION  TO. KUTTA  CONDITION 

A(KUTTA, J)  =  A ( KUTTA , J )  -  B 
A (KUTTA, KUTTA)  =  A ( KUTTA , KUTTA )  +  A(I,J) 

CONTINUE 


C 

C  FILL  IN  KNOWN  SIDES 

C 

A(I ,KUTTA+1 )  =  SINTHE (I )*COSALF  -  COSTHE ( I )*SINALF 
120  CONTINUE 

A ( KUTTA , KUTTA+ 1 )  =  -  ( COSTHE (1)  +  COSTHE (NODTOT) )*COSALF 
+  -  (SINTHE (1)  +  SINTHE (NODTOT) )*SINALF 

RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  SUBROUTINE  CORVOR  (SINALF,COSALF) 

C 

C  COMPUTE  THE  LOCAL  VELOCITIES  OF  CORE  VORTICES 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  CORVOR  (SINALF , COSALF) 

COMMON  /BOD/  IFLAG, NLOWER,NUPPER, NODTOT, X(202) ,Y(202) , 

+  COSTHE(201) ,SINTHE(201) ,SS ,NP1 ,NP2 

COMMON  'SING/  0(200) , GAMMA, OK (200) , GAMK 
COMMON  .  WAK /  TZ\l .  7OT ,  WAKE  ,  OT 

COMMON  , CORV/  CV( 200 ) , KC ( 200 ) , YC ( 200 ) ,M, ID , CC7X( 200 ) , CCVY(200 ) 

COMMON  /INF3/  AMY(200 , 201 ) , BMY(200 , 201 ) 

COMMON  /INF4/  CMY(200,200  ,CMX(200,200) 

COMMON  /POT/  PHI (200) ,PHIK( 200) 

COMMON  /GUST/  UG(200) .VG(200) ,XGF,UGUST,VGUST 
IF  (M  .EQ.  1)  GO  TO  46 
MM1  =  M  -  1 


C 

C  VELOCITY  COMPONENTS  OF  CORE  VORTICES  AT  CURRENT  TIME  STEP 

C 


UGC 

VGC 


=  0.0 
=  0.0 


111 


nnnnn 


DO  10 
XG 

IF  (XG 
UGC 
VGC 
VY 
+ 

vx 

+ 

DO  20 

VY 

VX 

CONTINUE 


N  =  1 ,MM1 

=  XC(N)*COSALF  +  YC(N)*SINALF 
.GT.  XGF)  GO  TO  5 
=  UGUST 
=  VGUST 

=  SS*BMY(N,NP1)*(GAMMA-GAMK)/WAXE+ 

( 1 . +UGC ) *S INALF+VGC*COSALF 
=  SS*AMY(N,NP1 )*(GAMMA-GAMK)/WAKE+ 

( 1 . +UGC ) ^COSALF- VGC*S INAL 
J  =  1 ,NODTOT 
=  VY  +  AMY(N,J)*QK(J)  +  BMY ( N , J ) *GAMK 
=  VX  -  BMY(N,J)*QK(J)  +  AMY(N, J) *GAMK 


ADD  CORE  VORTEX  CONTRIBUTION 

DO  30  MC  =  1 ,MM1 

IF  (MC  .EQ.  N)  GO  TO  30 

VY  =  VY  +  CMY(N,MC)*CV (MC) 

VX  =  VX  +  cmx(n,mc)*cv(mc) 

CONTINUE 

COORDINATES  OF  CORE  VORTICES  AT  NEXT  TIME  STEP 


(51 : 


10 

CONTINUE 

40 

CONTINUE 

50 

CONTINUE 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  SUBROUTINE  FANDM ( S INALF , COSALF )  C 

C  C 

C  COMPUTE  AND  PRINT  OUT  CD, CL, CM  C 

C  INTEGRATE  PRESSURE  DISTRIBUTION  BY  TRAPEZOIDAL  RULE  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  FANDM (S INALF, COSALF) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) , Y(202) , 

+  COSTHE(201) ,SINTHE(2Q1) ,SS,NP1,NP2 

COMMON  /CPD/  CP (200) 


CFX 

CFY 

CM 

DO  100 

XMID 

YMID 

DX 

DY 

CFX 

CFY 

CM 

CONTINUE 

CD 


=  0.0 
=  0.0 
=  0.0 

I  =  1 ,NODTOT 
=  .5*(x(I)  +  X(I+1)) 

=  .5*(Y  I)  +  Y  1+1)) 

=  X(I+1  -  XI 

=  Y(I+1)  -  Y  I 
=  CFX  +  CP ( I ) *DY 
=  CFY  -  CP(I)*DX 
=  CM  +  CP(I)*(DX*XMID  +  DY*YMID) 


=  CFX*COSALF  +  CFY*SINALF 
-  CFY*COSALF  -  CFX*S INALF 

( 6 , iooo )  :d,:l ,cm 


:M  =' ,710.3) 


WRITE  (6,1000)  ID, CL. CM 

1000  FORMAT ( / / , ’  CD  =',710.6,'  CL  =',710.5,'  CM  =',710. 5) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 

C  SUBROUTINE  GAUSS (NRHS)  C 

C  C 

C  SOLUTION  OF  LINEAR  ALGEBRAIC  SYSTEM  BY  C 

C  GAUSS  ELIMINATION  WITH  PARTIAL  PIVOTING  C 

C  C 

C  A  COEFFICIENT  MATRIX  C 

C  NEQNS  =  NUMBER  OF  EQUATIONS  C 


nnnnnnonnnmi  non  nnnn  non  nnrin  non  non non 


NRHS  =  NUMBER  OF  RIGHT  HAND  SIDES 


RIGHT-HAND  SIDES  AND  SOLUTIONS  STORED  IN 
COLUMNS  NEQNS+1  THRU  NEQNS +NRHS  OF  A 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  GAUSS (NRHS ) 

COMMON  /CDF /  A (20 1 , 211 ^ ,NEQNS 
NP  =  NEONS  -  1 

NTCT  =  NEONS  -  NRHS 

GAUSS  REDUCTION 

DO  150  I  =  2 ,NEQNS 

--  SEARCH  FOR  LARGEST  ENTRY  IN  (I-l)TH  COLUMN 
ON  OR  BELOW  MAIN  DIAGONAL 

IM  =1-1 

IMAX  =  IM 

AMAX  =  ABS (A(IM , IM) ) 

DO  110  J  =  I, NEONS 

IF  (AMAX  .GE.  A3S(A(J,IM)))  GO  TO  110 
IMAX  =  J 

AMAX  =  A3S (A( J . IM) ) 

110  CONTINUE 

SWITCH  (I-l)TH  AND  IMAXTH  EQUATIONS 

IF  (IMAX  .SO.  IM)  GO  TO  140 
00  120  J  =  IM.NTOT 

TEMP  =  A(IM,J; 

A  (,  IM ,  J )  =  A(IMAX.J) 

A(IMAX.J)  =  TEMP 
130  CONTINUE 

ELIMINATE  (I-l)TH  UNKNOWN  FROM 
ITH  THRU  (NEQNS)TH  EQUATIONS 

140  DO  150  J  =  I, NEQNS 

R  =  A( J, IM)/A(IM, IM) 

DO  150  K  =  I , NTOT 

150  A( J,K)  =  A( J,K)  -  R*A(IM,K) 

BACK  SUBSTITUTION 

DO  220  K  =  NP , NTOT 

A (NEQNS ,K)  =  A(NEQNS,K)/A(NEQNS, NEQNS) 

DO  210  L  =  2 , NEQNS 

I  =  NEQNS  +  1  -  L 

IP  =1+1 

J 


RETURN 


=  IP, NEQNS 

A(I,K  -  A(I,J)*A(J,K) 
A(I ,K)/A(I , I J 


DO  200 
200  A(I ,K) 
210  A(I ,K) 
220  CONTINUE 


SUBROUTINE  INDATA 

SET  PARAMETERS  OF  BODY  SHAPE 
FLOW  SITUATION,  AND  NODE  DISTRIBUTION 

USER  MUST  INPUT 

NLOWER  =  NUMBER  OF  NODES  ON  LOWER  SURFACE 
NUPPER  =  NUMBER  OF  NODES  ON  UPPER  SURFACE 
PLUS  DATA  ON  BODY  AND  SUBROUTINE  BODY 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


nnnnnnnnnn<<<>  n  no  n  o  n 


onnnnnn  onuonr 


SUBROUTINE  INDATA 
DIMENSION  TITLE (20) 

COMMON  /30D/  I  FLAG ,  NLOWER ,  NUPPER ,  NODTOT  , X ( 202 ) , Y ( 202 ) , 

+  COSTHS (201 ) , S INTHE ( 201 ) , SS ,NP1 ,NP2 

COMMON  /PAR/  NACA , TAU , ZPSMAX , PTMAX 
READ  (1.501)  ITITLE 
WRITS  (6,501)  ITITLE 
DO  10  I  =  1, ITITLE 
READ  .'1,502)  TITLE 

10  WRITE  (6,503)  TITLE 

501  FORMAT <31 5) 

502  FORMAT (20A4) 

503  FORMAT ( IX ,20A4) 

READ  (1,501)  I FLAG .NLOWER, NUPPER 

WRITE  (5.501)  IFLAG, NLOWER, NUPPER 

IF  (IFLAG  .NE.  0)  RETURN 

READ  (1,501)  NACA 

WRITE  (5,501)  NACA 

ISPS  =  NACA/lOuO 

I PTMAX  =  NACA/ 100  -  1Q*IEPS 

ITAU  =  NACA  -  10Q0*IEPS  -  100*IPTMAX 

EPSMAX  =  IEPS*C .01 

PTMAX  =  I PTMAX* 0 . 1 

TAU  =  ITAU*0 . 01 

IF  'ISPS  .LT.  10)  RETURN 

PTMAX  =  0.2025 

EPSMAX  =  2 . 6595*PTMAX**3 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  INFL  (MITR) 


CALCULATE  INFLUENCE  COEFFICIENTS 


ICCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  INFL  (NITR) 

COMMON  /BOD/  IFLAG .NLOWER .NUPPER .NODTOT , X( 202 ) ,Y(202) , 

+  COSTHE ( 201 ) , S INTHE (ZOl ) , SS ,NP1 ,NP2 

CCMMCN  /NUM/  PI , PI2INV 
COMMON  /WAX/  VYW , VXW . WAKE , DT 

COMMON  / CORV/  CV(200) ,XC(200) ,YC(200) ,M,TD,CCVX(200) ,CCVY(200) 
COMMON  /INF1/  AAN(201 , 201 ) ,BBN( 201 , 201 ) ,AYNP1 (201 ) , BYNP1 (201 j 
COMMON  /INF2/  CCN{201 ,200) ,CCT  201 ,200) ,CYNP1(200) ,CXNP1(200) 
COMMON  /INF3/  AMY(200 , 201 ) , BMY (200 , 201 ) 

COMMON  /INF4/  CMY(200 , 200 ) , CMX( 200 , 200 ) 

IF  ((M  .GT.  1)  .OR.  (NITR  .GT.  0))  GO  TO  510 

AAN(I, J)  :  NORMAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 

BY  UNIT  STRENGTH  DISTRIBUTED  SOURCE  ON  THE  J-TH  PANEL 

BBN(I , J)  :  NORMAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 

BY  UNIT  STRENGTH  DISTRIBUTED  VORTEX  ON  THE  J-TH  PANEL 


100 


DO  120  I  =  1, NODTOT 

XMID  =  .5*(xa)  +  X(I  +  1 ) ) 

V”’ID  -  .  z* ' " ‘  -  -  'T !  ~ ■' 

:p‘i:o  ;  =". .moStot 

FT  AN  =  P  V 

IF  (J  .EQ.  I)  GO  TO  100 
DXJ  =  XMID  -  X(J) 

DXJP  =  XMID  -  X  J+l) 

DYJ  =  YMID  -  Y(J) 

DYJP  =  YMID  -  Y( J+l) 

FLOG  =  . 5*ALOG( (DXJP*DXJP+DYJP*DYJP)/ (DXJ*DXJ+DYJ*DYJ) ) 
FTAN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE ( I ) *COSTHE ( J )  +  S INTHE ( I ) *S INTHE ( J ) 

STIMTJ  =  SINTHa ( I ) *COSTHE ( J)  -  COSTHE ( I ) *S INTHE ( J ) 
AAN(I,J)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ ) 


onnno 


110 

120 

510 


135 


C 

C 

C 

c 

c 

c 

c 

c 

c 


130 


140 


C 

C 

C 

C 

C 

C 

C 


BBN(I,J)  =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 

CONTINUE 
CONTINUE 
CONTINUE 
I  *  NP1 

XMID  =  . 5*( 

YMID  =  .  5*( 

DO  130  J  =  1, 

FLOG  =0.0 

FT AN  =  PI 

IF  (J  .EQ.  I)  GO  TO  135 
DXJ  =  XMID  -  X(J) 

DXJP  =  XMID  -  X  J+l) 

DYJ  =  YMID  -  Y(J) 

DYJP  =  YMID  -  Y( J+l) 

FLOG  =  .  5*ALOG(  (DXJP*DXJP+DYJP*DYJP)/ (DXJ*DXJ+DYJ*DYJ ) ) 

FTAN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE ( I ) *C0STHE ( J )  +  SINTHE(I)*SINTHE(J) 

STIMTJ  =  SINTHE ( I ) *COSTHE ( J)  -  COSTHE ( I ) *S INTHE ( J ) 

AAN ( I , J)  =  P I 2 INV* ( FTAN*CT IMT J  +  FLOG*STIMTJ) 

BBN ( I ,  J )  =  PI2INV* ( FLOG*CTIMT J  -  FTAN*STIMTJ) 

AYNPl(J)  s  Y  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
(NP1-TH  PANEL)  BY  UNIT  STRENGTH  DISTRIBUTED  SOURCE 
ON  J-TH  PANEL 

BYNPl(J)  s  Y  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
(NP1-TH  PANEL)  BY  UNIT  STRENGTH  DISTRIBUTED  VORTEX 
ON  J-TH  PANEL 


AYNPl(J)  =  ?I2INV*(FTAN*C0STHE(J)  -  r LOG*S INTHE ( J ) ' 

BYNP1 ( J)  =  ?I2INV*(FL0G*C0STHE(J)  +  FTAN*S INTHE ( J ) y 
CONTINUE 

DO  140  I  =  1 ,NODTOT 

XMID  =  .5*(X(I)  +  X(I+1)) 

YMID  =  .5*(Y(I)  +  Y(I+1) ) 

J  =  NP1 

DXJ  =  XMID  -  X(J) 

DXJP  =  XMID  -  X( J+l) 

DYJ  =  YMID  -  Y( J) 

DYJP  =  YMID  -  Y( J+l ) 

FLOG  =  . 5*ALOG((DXJP*DXJP+DYJP*DYJP)/ (DXJ*DXJ+DYJ*DYJ) ) 

FTAN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE (I)*COSTHE(J)  +  SINTHE (I )*SINTHE(J) 

STIMTJ  =  SINTHE ( I )*COSTHE(J)  -  COSTHE (I )*S INTHE (J) 

AAN(I.J)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 

BBN(I,J)  =  PI2INV*(FL0G*CT1MTJ  -  FTAN*STIMTJ) 

CONTINUE 

CCN ( I , J )  :  NORMAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 
BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 

CCT ( I , J )  :  TANGENTIAL  VELOCITY  INDUCED  AT  MID-POINT  OF  I-TH  PANEL 
BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 


IF  (M  .EQ.  1)  RETURN 

MM1  -  M  -  1 


.•  j.  . a  .  ji 


1)  SO  TO  520 
DO  220  I  =  l.NCDTOT 
XMID  =  O.S*(X(I)  +  X(I+1)) 

YMID  =  0.5*(Y(I)  +  Y  1+1)) 

DO  210  N  =  1.MM1 
DX  =  XMID  -  XC(N) 

DY  =  YMID  -  YC(N) 

DIST  =  SORT ( DX*DX+DY*DY ) 

COSTHN  =  DX/DIST 
SINTHN  =  DY/DIST 

CTIMTN  =  COSTHE (I )*COSTHN  +  SINTHE 
STIMTN  =  SINTHE ( I ) *COSTHN  -  COSTHE 
CCN(I,N)  *  -PI2INV*CTIMTN/DIST 


I)*SINTHN 
I )*SINTHN 


CCT(I,N)  =  -PI2INV*STIMTN/DIST 
210  CONTINUE 
220  CONTINUE 
520  CONTINUE 

I  =  NP1 

XMID  =  0 . 5*(X(I)  +  X(I+1)) 

YMID  =  0.5*(Y(I)  +  Y(I+1)) 

DO  230  N  =  1,MM1 
DX  =  XMID  -  XC(N) 

DY  =  YMID  -  YC(N) 

DIST  =  SORT ( DX*DX+DY*DY ) 

COSTHN  =  DX/DIST 
SINTKN  =  DY/DIST 

CTIMTN  =  COSTHE ( I ) *COSTHN-  +  SINTHE(I)*SINTHN 
STIMTN  =  S INTHE ( I ) *COSTHN  -  COSTHE ( I )*SINTHN 
CCN ( I , N )  =  -PI2INV*CTIMTN/DIST 
CCT(I ,N)  =  -PI2INV*STIMTN/DIST 
C 

C  CYNPl(N)  :  Y  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
C  (NP1-TH  PANEL)  BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 

C 

C  CXNPl(N)  ;  X  -  VELOCITY  INDUCED  AT  MID  POINT  OF  WAKE  ELEMENT 
C  (NP1-TH  PANEL)  BY  UNIT  STRENGTH  N-TH  CORE  VORTEX 

C 

CYNPl(N)  =  -PI2INV*COSTKN/DIST 
CXNPl(N)  =  +PI2INV*SINTHN/DIST 
230  CONTINUE 
C 

C  AMY(N, J)  :  Y  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
C  STRENGTH  DISTRIBUTED  SOURCE  ON  THE  J-TH  PANEL 


BMY(N, J) 


IF  (NITR 
DO  320  N 
XMID 
YMID 
DO  310 
DXJ 
DXJP 
DYJ 
DYJP 
FLOG 
FTAN 
AMY(N, J) 
BMY(N, J) 
CONTINUE 
CONTINUE 
CONTINUE 
DO  330 


Y  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  DISTRIBUTED  VORTEX  ON  THE  J-TH  PANEL 

.GT.  0)  GO  TO  530 
=  1 , MM1 
XC(N) 

YC(N) 

=  1 ,NODTOT 
XMID  -  X(J) 

XMID  -  X  J+l) 

YMID  -  Y(J) 

YMID  -  Y( J+l ) 

. 5  *ALOG ( ( DX JP*DX JP+DY JP*DY JP ) / ( DX J*DX J+DY J*DY J ) ) 
ATAN2 ( DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DY J ) 

=  PI2INV*(FTAN*COSTHE(J)  -  FLOG*SINTHE( J) ) 

=  PI2INV*(FLOG*COSTHE(J)  +  FTAN*SINTHE( J) ) 


FTAN*S INTHE ( 


XMID 

= 

XC(N) 

|  YMID 

= 

YC(N) 

J 

s 

NP1 

DXJ 

= 

XMID  - 

;  dxjp 

= 

XMID  - 

= 

YMID  - 

l  jYJ? 

= 

YMID  - 

r 

- 

.  5  rAL0G 

•(N, J)  = 
r(N,  J)  = 


=  ATAN2  (DYJP*DXJ-DXJP*DYJ ,  DXJP,'DXJ+DYJP’i;DYJ ) 
=  PI2INV*(FTAN*COSTHE ( J)  -  FLOG*SINTHE( J) ) 


CONTINUE 
CMY (N,MC) 


CMX(N,MC) 


PI2INV*(FLOG*COSTHE 


15); 


FTAN*S INTHE ( 


:  Y  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  MC-TH  CORE  VORTEX  OTHER  THAN  ITSELF 

:  X  -  VELOCITY  INDUCED  AT  N-TH  CORE  VORTEX  BY  UNIT 
STRENGTH  MC-TH  CORE  VORTEX  OTHER  THAN  ITSELF 
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no  on  ii  non  non  n  non  non 


410 

420 

CCt 


IF  (NITR  .GT.  0)  RETURN 


DO  420 

XMID 

YMID 

N  =  1.MM1 
=  KC ( N ) 

=  YC(N) 

DO  410 

MC  =  1 , MM1 

IF  (MC  , 

.EQ.  N)  GO  TO  410 

DX 

=  XMID  -  XC(MC) 

DY 

=  YMID  -  YC(MC) 

DIST 

=  SORT  ( DX*DX+DY’"DY ) 

CCSTHN 

=  OX/DIST 

SINTHN 

=  DY/DIST 

CMY(N,MC)  =  -  PI  2INV ’'CCSTHN /DIST 
CMX(N.MC)  =  +?I2INV*SINTHN/DIST 

CONTINUE 

CONTINUE 

RETURN 

END 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


SUBROUTINE  KUTTA  (ALPHA , SINALF , COS ALF , OMEGA, UX,UY) 


USING  KUTTA  CONDITION  TO  DETERMINE  VORTICITY 


ICCCC 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  KUTTA  ( ALPHA , 5 INALF . COSALF . OMEGA , UX , UY ) 

COMMON  /BOD/  IFLAG , MLCWER , NUPPER , NODTOT , K ( 202 ) ,Y(202) , 

+  CQSTHE (201 ) . SINTHE(2Q1) ,SS,NP1,NP2 

COMMON  /COF/  A( 201 ,211 ) ,NEQS 
COMMON  / SING/  0(200) , GAMMA (OK (200) ,GAMK 
COMMON  /WAK/  VYW , VXW . WAKE . DT 

COMMON  / CORV /  CV ( 200 ) , XC ( 200 ) , YC ( 200 ) , M , ID , CC7X (ZOOY,  C C7Y  ( 200 ) 
COMMON  / IMF 1 /  AAN( 201 , 201 ) , 33N( 201 , 201 ) ,AYNP1(201) ,3YNP1(201 ) 
COMMON  /INF2/  CCN!201 ,200) , CCT( 201 , 200 ) , CYNP1 (200) , CXNP1 (200 ) 
COMMON  /GUST/  UG ( 200  )  ,  7G (200),  XGF  ,  'JGUST ,  7GUST 
DIMENSION  81(200) ,32(200) ,AA(2) ,3B(2) 


RETRIEVE  SOLUTION  FROM  A-MATRIX 


50 


DO  50 
B1(I) 
B2  (I ) 


I  =  1, NODTOT 
=  A(I ,NP1 ) 

=  A(I ,NP2) 


FIND  VT  AT  TRAILING  EDGE  PANELS 


DO  130  K  *  1,2 

IF  (K  .EQ.  1)  I  =  1 

IF  (K  .EQ.  2)  I  =  NODTOT 

XMID  =  0.5  *  (X(I)  +  X(I+1)) 

YMID  =  0.5  *  (Y(I)  +  Y ( 1+1 ) ) 

VTANG  =  ( (1 .+UG(I))*COSALF-VG(I)*SINALF+UX)*COSTHE(I) 

+  +  ((l.+UG(I})*SINALF+VG(I)*COSALF+UY)*SINTHE(I) 

+  +  OMEGA* (YMID*COSTHE(I)  -  XMID*SINTHE (I ) ) 

AA(K)  =  -  AAN(I ,NP1)*SS/WAKE 

BB(K)  =  VTANG  +  AAN(I ,NP1 ) *SS*GAMMA/WAKE 

DO  120  J  =  1, NODTOT 

AA(K)  =  AA(X)  +  AANd,J)  -  BBN(I ,  J)*B1(J) 

3B(K)  =  3B(K)  -  3BN(I, J)*52(J) 

120  CONTINUE 

ADD  CORE  VORTEX  CONTRIBUTION 

IF  (M  .EQ.  1)  GO  TO  100 

MM1  =  M  -  1 

DO  110  N  =  1 ,MM1 

BB(K)  =  BB(K)  +  CCT(I,N)*CV(N) 

110  CONTINUE 
100  CONTINUE 
130  CONTINUE 

SATISFYING  KUTTA  CONDITION  --  SOLVE  FOR  VORTEX  STRENGTH 


117 


nnnnn 


c 


n 

s 


n 


* 


=  AA  ( 1 )  *AA  ( 1 )  -  AA  ('  2 '  *AA  ( 2 ) 

=  AAti  )*SB(  1 )  -  AA(2)’'BB{2)  -  SS/DT 
=  3B(  1  i*SB(  1 )  -  3B(2)*3BU)  +  2  .  *SS*GAMMA/DT 
=  3QRT(FF7'FF-EE*GG) 

=  (-FF  -  RADI) /EE 


C  CALCULATE  SOURCE  STRENGTH 

DO  150  I  =  1 , NODTOT 
160  QK ( I )  =  GAHKKB 1(1)  +  32(1) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  SUBROUTINE  NACA45 (2 , THICK, CAMSER , BETA)  C 

C  C 

C  EVALUATE  THICKNESS  AND  CAMBER  C 

C  FOR  NACA  4-  OR  5-DIGIT  AIRFOIL  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  NACA45 ( Z , THICK , CAMBER , BETA) 

COMMON  /PAR/  NACA , TAU , EPSMAX , PTMAX 
THICK  =0.0 

IF  iZ  .LT.  l.E-10)  GO  TO  100 

THICK  =  5 . *TAU*( . 2969^SQRT (2)  -  Z*(.126  +  Z*(.3537 
+  -  Z*( . 2343  -  Z*.1015)))) 

100  IF  (EPSMAX  .EQ.  0.0)  GO  TO  130 
IF  ( NACA  .GT.  9999)  GO  TO  140 
IF  (Z  .GT.  PTMAX)  GO  TO  110 

CAMBER  =  EPSMAX /PTMAX/ PTMAX*  ( 2 .  *?TMAX  -  Z-)*Z 
DCAMDX  =  2. *S?SMAX/ PTMAX/ PTMAX* (PTMAX  -  Z) 

GO  TO  120 

110  CAM3ER  =  EPSMAX/ (1. -PTMAX) **2*{1.  +  Z  -  2.*PTMAX)*(1.  -  Z) 

DCAMDX  =  2 .  ''EPSMAX/  ( 1 .  - PTMAX )**2* (PTMAX  -  Z) 

120  BETA  =  ATAN (DCAMDX) 

RETURN 

130  CAMBER  =0.0 

BETA  =0.0 

RETURN 

140  IF  (Z  .GT.  PTMAX)  GO  TO  150 

W  s  2 /PTMAX 

CAMBER  =  EPSMAX*W* ( (W  -  3.)*W  +  3.  -  PTMAX) 

DCAMDX  =  EPSMAX* 3.*W*(1.  -  W) /PTMAX 
GO  TO  120 

150  CAMBER  =  EPSMAX* (1.  -  Z) 

DCAMDX  =  -  EPSMAX 

GO  TO  120 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  SUBROUTINE  PRESS  ( S INALF , COSALF , OMEGA , UX , UY)  C 

C  C 

C  COMPUTE  UNSTEADY  FLOW  PRESSURE  DISTRIBUTION  C 

C  AND  VELOCITY  POTENTIAL  AT  MID-POINTS  OF  PANELS  C 

c  c 

-t  S*  ^  S*  ^  ^  r+  ~  r>  **  r*  r*  r*  r*  S'  r*  r*  **  r*  r*  f*  r*  C*  P  H  f*  r*  C*  C*  ^  ^ 

SUBROUTINE'*  press'"'’ SINALF^  COSALF  ,  OMEC-a"  Ux“jy’ 

COMMON  300 /  I FLAG , NLOWER . NUPPER , MODTCT , X ( 202 ) , Y( 202 ) , 
f  COSTHE(201),SINTHE(201),SS,NP1,NP2 

COMMON  / CPD/  CP (200) 

COMMON  /NUM/  PI,PI2INV 

COMMON  /SING/  0(200) , GAMMA, QK( 200) ,GAMK 

COMMON  /WAR/  VYW , VXW.WAKE , DT 

COMMON  / CORV/  CV(200) ,XC(200) ,YC(200) ,M,TD,CCVX(200) ,CCVY(200) 
COMMON  /INF1/  AAN(201 ,  201 ) , BBN( 201 , 201 ) , AYNP1 (201 )  , BYNP1 (201 ) 


COMMON  /INF1/  AAN(201 , 201 ) , BBN (201,201) , AYNP1 (201 ) , BYNP1 (201 ) 
COMMON  /INF2/  CCN(201 , 200) , CCT( 201 , 200 ) , CYNP1 (200) , CXNP1 (200 ) 
COMMON  /POT/  PHI(200),PHIK(200) 


COMMON  /GUST/  UG(200) ,VG(200) ,XGF,UGUST,VGUST 
COMMON  /EXTV/  UE(200) 


mRmmmmmm 


non  nnnnn  onn  non 


WRITE  (6,1000) 

FIND  TANGENTIAL  VELOCITY  VT  AT  MID-POINT  OF  I-TH  PANEL 


120 


DO  130 

XMID 

YMID 

DX 

DY 

DIST 

vsx 

VSY 

VS 

VTANG 


VTFREE 
VTANG 
DO  120 
VTANG 
CONTINUE 


I  =  1 ,NODTOT 
=  0.5  *  Ml)  +  X(I+1)) 

=  0.5  *  (Y(I  +  Y(I+1 ) ) 

:  %l\)  4 

=  SQRT(DX*DX+DY*DY) 

=  (1 . +UG(I ) )*COSALF-VG(I )*SINALF 
=  (l.+UG(I))*SINALF+VG  I  *C0SALF 
=  VSX*VSX  +  VSY*VSY 
=  ( ( 1 . +UG( I ) j *C0SALF-VG ( I ) *SINALF+UX ) *COSTHE (I) 

+  ( ( 1 . +UG ( I ) ) *  S INALF+VG ( I ) *C0  S ALF+UY ) *S I NTHE ( I ) 
+  OMEGA* (YMID*COSTHE ( I )  -  XMID*SINTHE ( I ) ) 

=  VTANG 

=  VTANG  +  SS*(GAMMA-GAMK)*AAN(I,NP1)/WAKE 
J  =  1 , N0DT0T 

=  VTANG  -  BBN(I, J)*QK(J)  +  AAN ( I , J ) *GAMK 


OKEGA*YMID  -  UX 
OMEGA*XMID  +  UY 


ADD  CORE  VORTEX  CONTRIBUTION 

IF  (M  .SO.  1)  GO  TO  150 
MM1  =  M  -  1 

DO  140  N  =  1,MM1 

VTANG  =  VTANG  +  CCT(I ,N)*CV(N) 

140  CONTINUE 
150  CONTINUE 

?HIK( I )  =  (VTANG- VTFREE ) *DIST 
CP (I)  =  VS  -  VTANG*VTANG 

UE ( I )  =  VTANG 

130  CONTINUE 


COMPUTE  DISTURBANCE  POTENTIAL  BY  LINE  INTEGRAL  OF  VELOCITY  FIELD 

INTEGRATION  FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 

NPHI  =  10  *  NLOWER 

PINK  =0.0 

XL  =  0.0 

DO  30  L  =  1 , NPHI 

FRACT  =  FLOAT (L)/FLOAT(NPHI) 

XLP  =  -10.0  *  (1.0  -  COS (0 . 5*PI*FRACT) ) 

DELX  =  XL  -  XLP 

XMID  =  0 . 5*(XL+XLP)*C05ALF 

YMID  =  0.5*(XL+XLP)*SINALF 

XL  =  XLP 

VELX  =  UGUST 

ADD  CONTRIBUTION  OF  J-TH  PANEL 

DO  20  J  =  1 ,NP1 

DXJ  =  XMID  -  X(J) 

DXJP  =  XMID  -  X(J+1) 

DYJ  =  YMID  -  Y(J) 

DYJ?  =  YMID  -  Y(J+1 : 

“LOG  =  .  5*AL0Gi  <  DXJ?*DXJ?-‘-DYJ?’'DYJ?  '  /  '  DXJ -DXJ  DY  J  *DY  I '  , 

FTAN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ ) 

CALMTJ  =  -COSALF*COSTHE ( J)  -  SINALF*SINTHE 7 J) 

SALMTJ  =  -SINALF*COSTHE ( J )  +  COSALF*SINTHE ( J) 

APY  =  PI2INV*(FTAN*CALMTJ  +  FLOG*SALMTJ) 

BPY  =  PI2INV* ( FLOG*CALMTJ  -  FTAN*SALMTJ ) 

IF  (J  .EQ.  NP1 )  GO  TO  40 
VELX  =  VELX  -  BPY*QK( J)  +GAMK*APY 

GO  TO  20 

40  VELX  =  VELX  +  S  S  *AP Y* ( GAMMA - GAMK ) / WAKE 
20  CONTINUE 
C 


C  ADD  CORE  VORTEX  CONTRIBUTION 

C 

IF  (M  .EQ.  1)  GO  TO  50 
MM1  =  M  -  1 
DO  60  N  =  1.MM1 
DX  =  XMID  -  XC(N) 

DY  =  YMID  -  YC(N) 

DIST  =  SORT ( DX*DX+DY*DY ) 

COSTHN  =  DX/DIST 
SIMTHN  =  DY/DIST 

SALMTN  =  -SINALF*COSTHN  +  COSALF^SINTHN 
CPT  =  -PI2INV*SALMTN/DIST 
60  VELX  =  VELX  +  CPT*CV(N) 

50  CONTINUE 

PINK  =  PINK  +  VELX  *  DELX 
30  CONTINUE 
C 

C  COMPUTE  DISTURBANCE  POTENTIAL  AT  MID-POINT  OF  I-TH  PANEL 

C 

C  LOWER  SURFACE 

C 

DO  230  I  =  1 , NLOWER 
PH  =  -PINK 

DO  240  J  =  I, NLOWER 
240  PH  =  PH  -  ?HIK( J) 

PHIK(I)  =  ?H 
230  CONTINUE 

DO  270  I  =  1, NLOWER- 1 

PHIK(I)  =  0.5*(PHIK(I)  +  PHIK(I+1) ) 

270  CONTINUE 

PHIK  ('NLOWER)  =  0 . 5*  (PHIK (NLOWER)  +  PINK) 

n 


C  UPPER  SURFACE 

c 

DO  250  I  =  NODTOT, NLOWER+1 , -1 
PH  =  -PINK 

DO  260  J  =  NLOWER+1 , I 
260  PH  *  PH  +  PHIK(J) 

PHIK ( I )  =  PH 
250  CONTINUE 

DO  280  I  =  NODTOT , NLOWER+2 , - 1 
PHIK(I)  =  0.5*(PHIK(I)  +  PHIK(I-l) ) 

280  CONTINUE 

PHIK (NLOWER+1)  =  0.5* (PHIK (NLOWER+1)  +  PINK) 
C 

C  COMPUTE  CP  AT  MID  POINT  OF  I-TH  PANEL 

C 


290 


DO  290  I  =  1, NODTOT 
XMID  =  . 5*(X(I )  +  X(I+1)) 

CP?I)  =  CP(I)  -  2.* (PHIK (I) -PHI (I) )/DT 
WRITE  (6,1050)  I ,XMID , QK( I ) , GAMK , CP ( I; 
CONTINUE 


UE(I) 


1000  FORMAT (/ ,4X, 1 J 1 ,4X, 1 X( J) 1 ,6X, 'Q(J) 1 ,5X, 'GAMMA' 


+  'CP(J) 1 , 6X, 'V 
1050  F0RMAT(I5 ,6F10 
RETURN 
END 


,/) 


5X, 


C  SUBROUTINE  SETUP 
C 

C  SETUP  COORDINATES  OF  PANEL  NODES  AND  SLOPES  OF  PANELS 

C  COORDINATES  ARE  READ  FROM  INPUT  DATA  FILE  UNLESS 

C  THE  AIRFOIL  IS  OF  NACA  XXXX  OR  NACA  230XX  TYPE 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  SETUP 

COMMON  / BOD/  IFLAG, NLOWER, NUPPER, NODTOT, X(202) ,Y(202) , 

+  C0STHE( 201 ) , SINTHE (201),SS,NP1,NP2 

COMMON  /NUM/  PI.PI2INV 


120 


nonnori. 


m 


m 


c 

c 


100 


110 


10 


501 

:o 


pi  = 

PI2INV  = 


3.1415926535 
■  5/PI 


SET  COORDINATES  OF  MODES  ON  30DY  SURFACE 

IF  ( I FLAG  .ME.  0)  GO  TO  10 

NPOINT  =  NLCWER 

SIGN  =  -1.0 

MSTART  =  0 

DO  11C  NSURF  =  1.2. 

DO  100  M  =  1 , NPOINT 

FRACT  =  FLOAT (N-i) /FLOAT (NPOINT) 

2  -  .5*(1.  -  COS (PI*FRACT ) ) 

I  =  NSTART  +  N 

CALL  30DY(Z , SIGN ,X( I ) ,Y( I ) ) 

CONTINUE 

NPOINT  =  NUPPER 
- x  GM  —  -  .  0 
NSTART  =  ML OWE R 
CONTINUE 

NQDTOT  =  NLOWER  +  NUPPER 
X(NODTOT+l )  =  X(l) 

Y(NODTOT+l)  =  Y { 1 ) 

GO  TO  20 

MCDTOT  =  MLOWER  -  NUPPER 
READ  (1,501  ( X  ( I ; ,1=1, MODTOT-1 ) 

WRITE  (6.501)  (X ( I ) , 1=1 , MODTOT+1 ) 

READ  (1,501:  (Y ( I) ,  1=1  .NODTOT-^1 ) 

WRITS  (6,501)  ,Y(I; , 1=1 ,NODTCT+l) 

FORMAT  ( 6F10 . 5 0 
MP1  =  MODTOT  -  1 


NP2 

NODTOT  -  2 

5 

t*  i 

1  SLOPES  OF 

PANELS  < 

SS 

0.0 

DO  200 

I 

=  1, NODTOT 

DX 

— 

X(I+1)  -  X! 

(I) 

DY 

= 

Y(I+1)  -  Yi 

I) 

DIST 

= 

SCRT(DX*DX 

+DY*DY) 

SS 

= 

S§  +  DIST 

SINTHE (1 

:) 

=  DY/DIST 

COSTHE ( 1 
CONTINUE 
RETURN 
END 

:) 

=  DX/DIST 

SUBROUTINE  TEWAK  (SINALF,COSALF) 

COMPUTE  WAKE  ELEMENT  AT  THE  TRAILING  EDGE 


200 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 
c 
c 
c 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  TEWAK  ( SINALF , COSALF) 

COMMON  /BOD/  IFLAG, NLOWER, NUPPER, NODTOT, X(202) ,Y(202) , 

+  COSTHE ( 201 ) , SINTHE ( 201 ) , SS , NP1 ,NP2 

COMMON  C OF/  A( 201 . 211 : .NEOS 
COMMON  SING/  C  :  ICC  >  ,  3AKMA  .  OK ( 200  >  .  C-AHK 
COMMON  WAX/  T/W , /MW . WAKE . DT 
COMMON  / WAK2/  7YWK , VXWK 

COMMON  / CORV/  CV( 200 ) , XC ( 200 ) , YC ( 200 ) , M ,TD , CCVX(200 ) , CCVY(200 ) 
COMMON  /IMF1/  AAN (201,201), BBN(201 , 201 ) , AYNP1 ( 201 ) , BYNP1 (201 ) 
COMMON  /INF2/  CCN( 201 , 200 ) , CCT( 201 , 200 ) , CYNF1 (200) , CXNP1 (200) 
COMMON  /GUST/  UG( 200 ) , VG(200 ) , XGF , UGUST , VGUST 
XMID  =  0.5  *  X(NP1)  +  X(NP2) ) 

YMID  =  0.5  *  (Y(NP1)  +  Y (NP2 ) ) 

UGW  =0.0 

VGW  =0.0 

XG  =  XMID*COSALF  +  YMID*SINALF 

IF  (XG  .GT.  XGF)  GO  TO  10 


121 


nonnn 


nnn  nnn  non  nnonnnnn  non 


UGW  =  UGUST 

VGW  =  '/GUST 

10  VYWK  =  (l.+UGW)*SINALF+VGW*COSALF 

VXWK  =  <"1.tUGW)*C0SALF-VGW*SINALF 

DO  120  J  =  1.N0DT0T 

VYWK  =  VYWK  +  AYN?1(J)*QK(J)  +  BYNP1 ( J)*GAMK 

120  VXWK  =  VXWK  -  BYNP1 ( J)XQK( J)  +  AYNP 1  ( J ) XGAMK 

ADD  CORE  VORTEX  CONTRIBUTION 

IF  (M  .20.  1)  GO  TO  140 

MM1  =  M  -  1 

DO  130  N  =  1 ,MM1 

VYWK  =  VYWK  +  CYNP1 (N)*CV(N) 

130  VXWK  =  VXWK  +  CXNP1(N)XCV(N) 

140  CONTINUE 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  VELDIS(SINALF,C0SALF) 

COMPUTE  STEADY  FLOW  PRESSURE  DISTRIBUTION 

AND  VELOCITY  POTENTIAL  AT  MID-POINTS  OF  PANELS 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  VELDIS (SINALF , COSALF) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPSR,NODTOT,X(202) ,Y(202) , 

+  COSTHE (201) , SINTHE (201 ) , SS ,NP1 ,NP2 

COMMON  /COF/  A(201 ,211)  , KUTTA 
COMMON  / CPD/  CP (200) 

COMMON  /NUM/  ?I,?I2INV 

COMMON  /SING/  0(200) .GAMMA, QK(200) ,GAMK 

COMMON  /'POT/’  PHI  (200  )  ,  ?HIK(  200 ) 

COMMON  /GUST/  UG( 200 ), VG( 200 ), XGF , UGUST ,VGUST 
COMMON  /EXTV/  UE(200) 

WRITE  (6,1000) 


RETRIEVE  SOLUTION  FROM  A-MATRIX 

DO  50  I  =  1 ,NODTOT 
50  Q(I)  =  A(I ,KUTTA+1) 

GAMMA  =  A ( KUTTA , KUTTA+1 ) 

FIND  VT  AND  CP  AT  MID-POINT  OF  I-TH  PANEL 


DO  130  I  =  1 ,NODTOT 

XMID  =  .5*(X(I)  +  X(I+1}) 

YMID  =  . 5*(Y( I )  +  Y ( 1+1 ) ) 

VTANG  =  COS ALF*COSTHE ( I )  +  SINALF*SINTHE(I) 
VTFREE  =  VTANG 


ADD  CONTRIBUTION  OF  J-TH  PANEL 


i 

! 

! 


100 


120 

i 


DO  120  J  =  1 , NODTOT 

FLOG  =0.0 

FTAN  =  PI 

_?  . J  . EO .  . i  GO  TO  100 
DXJ  =  :1MID  -  X(J) 

DXJP  =  XMID  -  X(J+1) 

DYJ  =  YMID  -  Y(J) 

DYJP  =  YMID  -  Y(J+1) 

FLOG  =  .5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

FTAN  =  ATAN2 (DYJP*DXJ-DXJP*DYJ ,DXJPxDXJ+DYJP*DYJ ) 

CTIMTJ  =  COSTHE ( I )*COSTHE(J)  +  SINTHE(I)*SINTHE(J) 

STIMTJ  =  SINTHE (I )*COSTHE(J)  -  COSTHE (I )*SINTHE( J) 

AA  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 

B  =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 

VTANG  =  VTANG  -  B^Q^J)  +GAMMA*AA 

CONTINUE 


nnnnnn 


non  onnnn  nnn  ooonn  non 


130 


20 

30 


CP(I) 

UE(I) 


=  1.0  -  VTANG*VTANG 
=  VTANG 

WRITS  (6,1050)  I , XMID ,  Q  ( I )  , GAMMA , CP ( I ) , UE ( I ) 


INITIAL  SET-UP  FOR  DISTURBANCE  POTENTIAL  CALCULATION 


^  :  fiij 

=  SQRT(DX*DX+DY*DY) 

-  /TrTla^tr!_^rpcDE■c•  \  *nT< 


DX 
DY 

DIST 

PHI (I)  -  (0TANG-VTFREE)*DIST 

CONTINUE 


COMPUTE  DISTURBANCE  POTENTIAL  BY  LINE  INTEGRAL  OF  VELOCITY  FIELD 
INTEGRATION  FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 


NPHI 

PIN 

XL 

DO  30 

FRACT 

XLP 

DELX 

XMID 

YMID 

XL 

VELX 


=  10  *  NLOWER 
=  0.0 
=  0.0 

L  =  1 , NPHI 

=  FLOAT (L) /FLOAT (NPHI) 

=  -10.0  *  (1.0  -  COS (6 . 5*PI*FRACT) ) 
s  XL  -  XLP 

=  0 . 5*(XL+XLP)*C0SALF 
=  0 . 5*(XL+XLP)*SINALF 
=  XLP 
=  UGUST 

ADD  CONTRIBUTION  OF  J-TH  PANEL 


DO  20 

DXJ 

DXJP 

DYJ 

DYJP 

FLOG 

FTAN 

CALMTJ 

SALMTJ 

APY 

BPY 

VELX 

CONTINUE 

PIN 

CONTINUE 


J  =  1 , NODTOT 
=  XMID  -  X(J) 

=  XMID  -  X(J+1) 
=  YMID  -  Y(J) 

=  YMID  -  Y  J+l) 


=  . 5  *ALOG ( ( DXJP*DX JP+DY JP*DYJP ) / (DX J*DX J+DY J*DY J ) ) 
=  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

=  -COSALF*COSTHE ( J )  -  SINALF*SINTHE(J) 

=  -SINALF*COSTHE(J)  +  COSALF*SINTHE( J) 

=  PI2INV*(FTAN*CALMTJ  +  FLOG*SALMTJ) 

=  PI2INV*(FLOG* C ALMT J  -  FTAN*SALMTJ) 

=  VELX  -  BPY*Q( J)  +GAMMA*APY 


=  PIN  +  VELX  *  DELX 


240 

230 


COMPUTE  DISTURBANCE  POTENTIAL  AT  MID-POINT  OF  I-TH  PANEL 

LOWER  SURFACE 

DO  230  I  =  1 , NLOWER 
PH  =  -PIN 

DO  240  J  =  I, NLOWER 
PH  =  PH  -  PHI(J) 

PHI (I)  =  PH 

CONTINUE 

DO  270  I  =  1. NLOWER- 1 

PHI  (I)  =  0.5*'  PHI;  I)  -r  ?HI(I+I)) 

CONTINUE 

PHI (NLOWER)  =  0.5* (PHI (NLOWER)  +  PIN) 

UPPER  SURFACE 

DO  250  I  =  NODTOT, NLOWER+1,-1 
PH  =  -PIN 

DO  260  J  =  NLOWER+1 , I 
PH  *  PH  +  PHI(J) 

PHI (I)  =  PH 

CONTINUE 

DO  280  I  =  NODTOT, NLOWER+2,-1 


230 

1000 

1050 


PHI(I)  =  0.5*(PHI(I)  +  PHI(I-l))  ' 

CONTINUE 

PHI (NLOWER+1)  =  0 . 5* (PHI (NLOWER+1 )  +  PIN) 

FORMAT ( 7 , 4X , ‘ J 1 ,4X, 'X(J) 1 ,6X, 'Q(J) ' ,5X, ‘GAMMA' ,5X, 
+  1 CP(J) 1 ,6X, 'V(j)‘ ,/) 

FORMAT ( 15 , 5F10 . b ) 


RETURN 
END 


124 


APPENDIX  B 

EXAMPLE  INPUT  DATA  FOR  PROGRAM  U2DIIF 


THIS  IS  AN  EXAMPLE  OUTPUT  DATA  OBTAINABLE  FROM  PROGRAM  U2DIIF 
AIRFOIL  :  MISE5  8.4%  THICKNESS  (COORDINATES  ARE  INPUT  BY  USER) 

PANEL  NUMBER  :  NLOWER  =  25  ,  NUPPER  =25 

AIRFOIL  MOTION  :  MODIFIED  RAMP  AOA  CHANGE  ABOUT  MID  CHORD 

INITIAL  AOA  :  2.5  DEGREES 

FINAL  AOA  :  7.5  DEGREES 

AOA  RISE  TIME  :  1.5  CHORD  LENGTHS 

COMPUTATION  TIME  STEP  :  0.05  DURING  TRANSIENT  MOTION,  INCREASES 

PROGRESSIVELY  AFTER  FINAL  AOA  IS  REACHED 
*********^ *********2 ******* **3 ******** *4* ******** 5 *********6 *********7 


01  25 

25 

1.000000 

0.994853 

0.980866 

0.958884 

0.929536 

0.893455 

0.351308 

0.803815 

0. ''51753 

0.695943 

0.637271 

0.576620 

0.514913 

0.453098 

0.392084 

0.332794 

0.276105 

0.222865 

0.173361 

0.129819 

0.091393 

0.059146 

0.033560 

0.015010 

0.0C3767 

0.000000 

0.003767 

0.015003 

0.033560 

0.059146 

0.091393 

0.129319 

0.173861 

0.222865 

0.276105 

0.332791 

0.392032 

0.453G95 

0.514915 

0.576617 

0.637266 

0.695946 

0.751750 

0.803815 

0.851308 

0.393455 

0.929536 

0.958884 

0.980366 

0.994858 

1.000000 

0.000000 

-0.000782 

-0.002784 

-0.005721 

-0.009351 

-0.013459 

-0.017837 

-0.022235 

-0.026613 

-0.030671 

-0.034289 

-0.037341 

-0.039712 

-0.041314 

-0.042083 

-0.041979 

-0.040979 

-0.039096 

-0.036360 

-0.032820 

-0 . 028555 

-0.023651 

-0.013220 

-0.012379 

-0.006259 

0.000000 

0.006259 

0.012379 

0.018220 

0.023651 

0.028555 

0.032820 

0.036360 

0.039096 

0.040979 

0.041979 

0.042083 

0.041314 

0.039712 

0.037341 

0.034289 

0.030671 

0.026618 

0.022285 

0.017837 

0.013459 

0.009351 

0.005721 

0.002784 

0.000782 

0.000000 

2.50000 

5.000000 

1.5 

0.0 

0.5 

0.0 

0.000000 

0.000000 

0.000000 

2.000000 

0.050000 

0.0001 

0.000000 

APPENDIX  C 

EXAMPLE  OUTPUT  DATA  FROM  PROGRAM  U2DIIF 


DATA  READ  FROM  FILE  CODE  1 


11 

***x**rt**^**Ti(K**'j>rxx2*'*:*******3*********4***51t*****5*********5******,,t**7 
THIS  IS  AN  EXAMPLE  OUTPUT  DATA  OBTAINABLE  FROM  PROGRAM  U2DIIF 
AIRFOIL  :  MISES  8.4%  THICKNESS  (COORDINATES  ARE  INPUT  BY  USER) 

PANEL  NUMBER  :  NLOWER  =  25  ,  NUPPER  =25 

AIRFOIL  MOTION  :  MODIFIED  RAMP  AOA  CHANGE  ABOUT  MID  CHORD 

INITIAL  AOA  j  2.5  DEGREES 

FINAL  AOA  :  7.5  DEGREES 

ACA  RISE  TIME  :  1.5  CHORD  LENGTHS 

COMPUTATION  TIME  STEP  :  0.05  DURING  TRANSIENT  MOTION,  INCREASES 

PROGRESSIVELY  AFTER  FINAL  AOA  IS  REACHED 

1  25  25 

1.000000  0.994353  0.980366  0.958884  0.929536  0.893455 

0.351308  0.803815  0.751753  0.695948  0.637271  0.576620 

0.514918  0.453098  0.392084  0.332794  0.276105  0.222865 

0.173361  0.129819  0.091393  0.059146  0.033560  0.015010 

0.003767  0.000000  0.003767  0.015008  0.033560  0.059146 

0.091393  0.129819  0.173861  0.222865  0.276105  0.332791 

0.392082  0.453095  0.514915  0.576617  0.637266  0.695946 

0.751750  0.303815  0.851308  0.893455  0.929536  0.958884 

0.980866  0.994858  1.000000 

0.000000  -0.000782  -0.002784  -0.005721  -0.009351  -0.013459 

-0.017837  -0.022285  -0.026618  -0.030671  -0.034289  -0.037341 

-0.039712  -0.041314  -0.042083  -0.041979  -0.040979  -0.039096 

-0.036360  -0.032320  -0.028555  -0.023651  -0.018220  -0.012379 

-0.006259  0.000000  0.006259  0.012379  0.018220  0.023651 

0.028555  0.032820  0.036360  0.039096  0.040979  0.041979 

0.042083  0.041314  0.039712  0.037341  0.034289  0.030671 

0.026618  0.022285  0.017837  0.013459  0.009351  0.005721 

0.002784  0.000782  0.000000 

2.500000  5.000000  1.500000  0.000000  0.500000  0.000000  0.000000 

0.000000  0.000000  0.000000 

2.000000  0.050000  0.000100  0.000000 


AIRFOIL  PERIMETER  LENGTH  =  2.018599 


STEADY  FLOW  SOLUTION  AT  ALPHA  =  2.500000 


T 

X(J) 

Q(  J) 

i 

0.997429 

0.355723 

2 

0.987862 

0.356105 

3 

0.969875 

0.365026 

4 

0.944210 

0.378836 

5 

0.911495 

0.394973 

6 

0.872381 

0.412926 

7 

0.827561 

0.432568 

8 

0.777784 

0.453710 

9 

0.723850 

0.476112 

10 

0.666609 

0.500047 

11 

0.606945 

0.525455 

12 

0.545769 

0.552654 

GAMMA  C?(J)  7(J) 

0.074003  0.316305  -0.326359 
0.074003  0.206074  -0.891025 
0.074003  0.133790  -0.930704 
0.074003  0.082276  -0.957979 
0.074003  0.043033  -0.978247 
0.074003  0.012034  -0.993965 
0.074003  -0.012724  -1.006342 
0.074003  -0.032414  -1.016078 
0.074003  -0.048010  -1.023724 
0.074003  -0.059905  -1.029517 
0.074003  -0.068405  -1.033637 
0.074003  -0.073563  -1.036129 


13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 


0.484008 

0.422591 

0.362439 

0.304449 

0.249485 

0.198363 

0.151840 

0.110606 

0.075269 

0.046353 

0.024285 

0.009388 

0.001884 

0.001884 

0.009387 

0.024234 

0.046353 

0.075269 

0.110606 

0.151840 

0.198363 

0.249485 

0.304448 

0.362436 

0.422588 

0.484005 

0.545766 

0.606941 

0.666606 

0.723848 

0.777732 

0.827561 

0.872381 

0.911495 

0.944210 

0.969875 

0.987862 

0.997429 


0.581715 

0.612882 

0.646430 

0.683297 

0.723595 

0.768523 

0.819732 

0.879132 

0.951277 

1.043426 

1.172784 

1.380641 

1.653644 

0.545367 

-0.275210 

-0.497235 

-0.580394 

-0.618032 

-0.636699 

-0.645594 

-0.649755 

-0.650971 

-0.650596 

-0.649624 

-0.648020 

-0.646281 

-0.644572 

-0.642990 

-0.641396 

-0.639980 

-0.638504 

-0.637214 

-0.635880 

-0.634260 

-0.632196 

-0.629793 

-0.625960 

-0.619834 


0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 

0.074003 


-0.075309 

-0.073531 

-0.067928 

-0.057668 

-0.041891 

-0.019114 

0.013374 

0.059669 

0.127774 

0.232984 

0.409236 

0.727179 

0.945112 

-0.815326 

-1.084184 

-0.872601 

-0.723499 

-0.624167 

-0.552954 

-0.498371 

-0.453794 

-0.415592 

-0.381557 

-0.349957 

-0.319875 

-0.290579 

-0.261352 

-0.231604 

-0.200941 

-0.168763 

-0.134640 

-0.097676 

-0.056864 

-0.010900 

0.042413 

0.107609 

0.193062 

0.316307 


-1.036971 

-1.036114 

-1.033406 

-1.028430 

-1.020731 

-1.009512 

-0.993290 

-0.969707 

-0.933931 

-0.375795 

-0.763612 

-0.522323 

0.234282 

•1.347341 

1.443670 

1.368430 

1.312821 

1.274428 

1.246176 

1.224080 

1.205734 

1.189787 

1.175397 

1.161877 

1.148858 

1.136037 

1.123099 

1.109776 

1.095875 

1.081094 

1.065195 

1.047701 

1.028039 

1.005436 

0.978564 

0.944665 

0.898297 

0.826858 


CD  *  0.000829  CL  *  0.303076  CM  =  -0.080325 


*************************************** 

***  BEGIN  UNSTEADY  FLOW  SOLUTION  **** 
*************************************** 


TIME  STEP  TK  =  0.050000 


TK  -  TKM1  *  0.050000 


ALPHA (T)  =  2.516295  OMEGA (T)  =  -0.011248 

U(T)  =  0.000000  V(T)  =  -0.005624 


NITR  VXW  VYW 

0  0.999048  0.043619 

1  0.907832  0.005991 

2  0.904138  0.007297 

3  0.903985  0.007241 


WAKE  THETA 

0.050000  0.043633 
0.045393  0.006600 
0.045208  0.008070 
0.045201  0.008010 


GAMMA 

0.740032E-01 
0.744799E-01 
0.744662E-01 
0 . 744652E-01 


CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  ■  3 

J  X(J)  Q(J)  GAMMA  CP(J)  V(J) 

1  0.997429  0.435837  0.074466  0.299470  -0.838202 


2 

0.987862 

0.422852 

0.074466 

0.196120 

-0.399569 

3 

0.969875 

0.421798 

0.074466 

0.132077 

-0.936976 

4 

0.944210 

0.427256 

0.074466 

0.083439 

-0.962416 

5 

0.911495 

0.435974 

0.074466 

0.055859 

-0.931131 

6 

0.872381 

0.447131 

0.074466 

0.029322 

-0.995676 

7 

0.827561 

0.460598 

0.074466 

0.008149 

-1.007027 

8 

0.777784 

0.475913 

0.074466  -0.010442 

-1.015962 

9 

0.723850 

0.492830 

0.074466  -0.026773 

-1.022942 

10 

0.666609 

0.511571 

0.074466  -0.041132 

-1.023223 

11 

0.606945 

0.532044 

0.074466  -0.053494 

-1.032001 

12 

0.545769 

0.554564 

0.074466  -0.063531 

-1.034262 

13 

0.484008 

0.579201 

0.074466  -0.070990 

-1.035013 

14 

0.422591 

0.606206 

0.074466  -0.075185 

-1.034204 

15 

0.362439 

0.635915 

0.074466  -0.075487 

-1.031672 

16 

0.304449 

0.669133 

0.074466  -0.070722 

-1.027008 

17 

0.249485 

0.706140 

0.074466  -0.059745 

-1.019772 

18 

0.198363 

0.748100 

0.074466  -0.040826 

-1.009171 

19 

0.151840 

0.796683 

0.074466  -0.011154 

-0.992762 

20 

0.110606 

0.853825 

0.074466 

0.033396 

-0.971231 

21 

0.075269 

0.924099 

0.074466 

0.100715 

-0.936353 

22 

0.046353 

1.014820 

0.074466 

0.205876 

-0.880676 

23 

0.024285 

1.143358 

0.074466 

0.382654 

-0.776542 

24 

0.009388 

1.351849 

0.074466 

0.703606 

-0.535371 

25 

0.001884 

1.634712 

0.074466 

0.952740 

0.209596 

26 

0.001834 

0.564272 

0.074466  -0.746972 

1.322541 

27 

0.009387  -0.246433 

0.074466  -1.037466 

1.430099 

28 

0.024284  -0.467818 

0.074466  -0.838061 

1.360476 

29 

0.046353  -0.551795 

0.074466  -0.693563 

1.307915 

30 

0.075269  -0.590860 

0.074466  -0.596473 

1.271481 

31 

0.110606  -0.611392 

0.074466  -0.527117 

1.244626 

32 

0.151840  -0.622549 

0.074466  -0.474313 

1.223530 

33 

0.198363  -0.629333 

0.074466  -0.433321 

1 .206047 

34 

0.249485  -0.633515 

0.074466  -0.399064 

1.190716 

35 

0.304448  -0.636433 

0.074466  -0.369826 

1 . 176793 

36 

0.362436  -0.639059 

0.074466  -0.343641 

1.163536 

37 

0.422588  -0.641343 

0.074466  -0.319346 

1.150745 

38 

0.484005  -0.643768 

0.074466  -0.295850 

1.137963 

39 

0.545766  -0.646483 

0.074466  -0.272088 

1.124936 

40 

0.606941  -0.649578 

0.074466  -0.247068 

1.111386 

41 

0.666606  -0.652918 

0.074466  -0.220048 

1.097122 

42 

0.723848  -0.656699 

0.074466  -0.190134 

1.081849 

43 

0.777782  -0.660707 

0.074466  -0.156552 

1.065285 

44 

0.827561  -0.665236 

0.074466  -0.118313 

1.046980 

45 

0.872381  -0.670135 

0.074466  -0.074279 

1.026307 

46 

0.911495  -0.675261 

0.074466  -0.023233 

1.002476 

47 

0.944210  -0.680614 

0.074466 

0.036802 

0.974108 

43 

0.969875  -0.686543 

0.074466 

0.109850 

0.938386 

49 

0.987862  -0.692719 

0.074466 

0.203356 

0.889792 

50 

0.997429  -0.699982 

0.074466 

0.333160 

0.815602 

CD  ■  0.001539  CL  =  0.302054  CM  =  -0.088450 


TRAILING  VORTICES  DATA 

M  X(M)  Y(M)  CIRC 

1  1.022599  0.9C0131  -0.300933 


TIME  STEP  TK  =  0.749999  TK  -  TKM1  =  0.050000 

ALPHA (T)  =  4.999996  OMEGA (T)  =  -0.087266 

U(T)  =  0.000000  V(T)  =  -0.043633 

NITR  VXW  VYW  WAKE  THETA  GAMMA 

0  0.905684  -0.000916  0.045284  -0.001012  0.103235E+00 


1  0.905735  -0.000649  0.045287  -0.000717 


0 . 106563E+00 


CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =  1 


J 

X(J)  2<J5 

GAMMA 

CP(J) 

0.997429  1.115997 

0.106565 

0.311649 

• 

2 

0.987862  1.060332 

0.106565 

0.221106 

- 

2 

0.969373  1.031653 

0 . 106555 

0.170306 

- 

4 

0.944210  1.012085 

0 .106565 

0.:41153 

- 

5 

0.911495  0.993141 

0.106565 

0.123949 

- 

6 

0.372331  0.985265 

0 . 106565 

0.114073 

- 

7 

0.327561  0.974241 

0.106565 

0.109016 

- 

8 

0.777784  0.964993 

0.106565 

0.107136 

- 

9 

0.723850  0.957451 

0.106565 

0.107130 

- 

10 

0.566609  0.952120 

0.106565 

0.108453 

- 

*  X 

0.506945  0.949235 

0 . 106565 

0.110657 

- 

12 

0.545^69  0.949293 

0.106565 

0.113665 

- 

i  t 

0.484008  0.952960 

0.106565 

0.117540 

- 

14 

0.422591  0.960451 

0.106565 

0.122430 

- 

15 

0.362439  0.972455 

0.106565 

0.128967 

- 

16 

0.304449  0.990059 

0.106565 

0.138078 

-1 

17 

0.249485  1.013807 

0.106565 

0.150962 

13 

0.198263  1.045081 

0.106565 

0.169616 

- 

19 

0.151340  1.035735 

0 . 106565 

0.197364 

- 

20 

0.110606  1.123024 

0 . 106565 

0.239070 

— 

21 

0.075269  1.206453 

0 . 106565 

0.303822 

- 

22 

0.046253  1.298127 

0.106565 

0.407818 

- 

23 

0.024285  1.423849 

0 . 106565 

0.533190 

- 

24 

0.009333  1.531804 

0.106565 

0.872165 

- 

25 

0.001334  1.319807 

0.106565 

0.765241 

26 

0.001834  0.3731^9 

0.106565 

-1.559307 

0.009237  -0.529374 

0.106565 

-1 . 564167 

23 

0.124234  -0.755145 

0.106565 

-1.202155 

29 

0.046353  -0.326350 

0 .106565 

-0.991312 

30 

0.075269  -0.374103 

0.106565 

-0.864650 

31 

0.110606  -0.896239 

0.106565 

-0.781035 

32 

0.151340  -0.912C96 

0.106565 

-0.721021 

33 

0.198363  -0.926607 

0.106565 

-0.674005 

34 

0.249435  -0.941343 

0.106565 

-0.634257 

35 

0.304448  -0.957406 

0 . 106565 

-0.598321 

35 

0.362436  -0.975544 

0.106565 

-0.563539 

37 

0.422588  -0.995441 

0.106565 

-0.528592 

38 

0.484005  -1.017302 

0.106565 

-0.492391 

39 

0.545766  -1.041010 

0.106565 

-0.454043 

40 

0.606941  -1.066387 

0.106565 

-0.412924 

41 

0.666606  -1.093023 

0.106565 

-0.368740 

42 

0.723848  -1.120818 

0.106565 

-0.321026 

43 

0.777782  -1.149241 

0.106565 

-0.269645 

44 

0.827561  -1.178310 

0.106565 

-0.213999 

45 

0.872381  -1.207644 

0.106565 

-0.153563 

46 

0.911495  -1.236895 

0.106565 

-0.087684 

47 

0.944210  -1.265976 

0.106565 

-0.014749 

48 

0.969875  -1.296105 

0.106565 

0.069103 

49 

0.987862  -1.330154 

0.106565 

0.171240 

50 

0.997429  -1.380203 

0.106565 

0.308161 

:d 

=  1. 333337  :L  : 

-  O.o4533i 

3  CM  = 

-0 

TRAILING  VORTICES  DATA 

M 

X(M)  Y<M) 

CIRC 

1 

1.700760  0.077052 

-0.C00933 

2 

1.651289  0.072156 

-0.001625 

3 

1.602002  0.066331 

-0.002229 

4 

1.552845  0.060031 

-0.002784 

5 

1.503779  0.053471 

-0.003304 

6 

1.454797  0.046788 

-0.003797 

7 

1.405895  0.040107 

-0.004256 

V(J) 

0.908159 

0.957033 

0.983709 

0.993976 

1.0C8098 

1.013416 

1.016142 

1.017013 

1.016594 

1.015089 

1.012569 

1.009317 

1.004971 

0.999505 

0.992654 

0.983855 

0.972486 

0.957451 

0.936848 

0.907675 

0.363894 

0.793054 

0.663132 

0.369686 

0.485826 

1.595323 

1.590937 

1.468068 

1.389539 

1.338450 

1.302194 

1.274529 

1.251843 

1.232132 

1.214148 

1.196886 

1.179831 

1.162518 

1.144540 

1.125555 

1.105335 

1.083543 

1.059939 

1.034027 

1.005269 

0.972965 

0.935775 

0.890819 

0.832327 

0.746084 


.224293 


129 


3 

1.357057 

0.033554  -0.004637 

a 

1.303305 

0.027216  -0.005100 

10 

1.259669 

0.021159  -0.005481 

•  i 

1.211173 

0.015541  -0.005301 

■  i 

1.162916 

0.010431  -0.006100 

13 

1.115043 

0.006079  -0.006349 

14 

1.067913 

0.002406  -0.006559 

15 

1.022643 

-0 .000016  -0.006720 

TIME  STEP  TK  =  1.449992 


TK  -  TKM1  =  0.050000 


AL?HA(7)  =  7.433698  OMEGA(T)  =  -0.011249 

U(T)  =  0.000000  V(T)  =  -0.005625 


NITR  7X W 


7YW  WAKE  THETA  GAMMA 


0  C. 901699  0.009872 

1  0.901200  0.011028 


0.045088  0.010948  0.145377E+00 

0.045063  0.012236  0.146997E+00 


CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =  1 


J 

Vf 

V  -  / 

QU) 

GAMMA 

C?(J) 

V(J) 

i 

0.997429 

1.101898 

0.146996 

0.332408  -0.864078 

0.937362 

1.099609 

0.146996 

0.228004  -0.921188 

3 

0.969875 

1.116333 

0.146996 

0.160643  -0.955077 

4 

0.944210 

1.140985 

0.146996 

0.114868  -0.9765.79 

■> 

0.911495 

1 .  .53619 

0.146996 

0.032627  -0.990894 

5 

0.372331 

1.197904 

0.146996 

0.060528  -1.000257 

■7 

0  .  SC'Sol 

1.228671 

0.146996 

0.046710  -1.005878 

a 

0.777734 

1.260826 

0.146996 

0.039936  -1.008490 

9 

0.723850 

1.294203 

0.146996 

0.039112  -1.008660 

10 

0.666609 

1.329222 

0.146996 

0.043766  -1.006560 

11 

0.606945 

1.366027 

0.146996 

0.053332  -1.002352 

12 

0.545759 

1.405159 

0.146996 

0.067653  -0.995947 

13 

0.484008 

1.446882 

0.146996 

0.086578  -0.987214 

14 

0.422591 

1.491619 

0.146996 

0.110103  -0.975912 

15 

0.362439 

1.539856 

0.146996 

0.138557  -0.961606 

16 

0.304449 

1.592619 

0.146996 

0.172964  -0.943442 

17 

0.249485 

1.650390 

0.146996 

0.214399  -0.920465 

18 

0.198363 

1.714437 

0.146996 

0.265094  -0.890942 

19 

0.151840 

1.786607 

0.146996 

0.328746  -0.851951 

20 

0.110606 

1.868882 

0.146996 

0.410545  -0.798852 

21 

0.075269 

1.965495 

0.146996 

0.519597  -0.722328 

22 

0.046353 

2.082446 

0.146996 

0.668455  -0.603478 

23 

0.024285 

2.230950 

0.146996 

0.866926  -0.394782 

24 

0.009388 

2.419835 

0.146996 

1.009466 

0.050859 

25 

0.001884 

2.339782 

0.146996 

-0.471799 

1.214887 

26 

0.001834 

-0.156346 

0.146996 

-4.389847 

2.320095 

27 

0.009387 

-1.322127 

0.146996 

-3.033496 

2.003043 

23 

0.024284 

-1.560176 

0.146996 

-2.015200 

1.727204 

29 

0.046353 

-1.622657 

0.146996 

-1.505832 

1.569752 

■0 

j . 075259 

-1.634560 

0.146996 

-1.212790 

1.470549 

3 1 

j . .13306 

-1 . 523102 

0 . .46996 

-1.021304 

1.401554 

2  Z 

j. .31340 

- ; .  112625 

j . 146996 

-0.335616 

1.350007 

j  2 

3 . 1 93253 

-..596423 

J. -46996 

-0 .730655 

1.309004 

34 

0.249485 

-1.578180 

0.146996 

-0.695048 

1.274895 

35 

0.304448 

-1.560042 

0.146996 

-0.621833 

1.245407 

36 

0.362436 

-1.542861 

0.146996 

-0.556431 

1.218916 

37 

0.422588 

-1.526391 

0.146996 

-0.496616 

1.194569 

38 

0.484005 

-1.510876 

0.146996 

-0.440592 

1.171595 

39 

0.545766 

-1.496317 

0.146996 

-0.387199 

1.149437 

40 

0.606941 

-1.482638 

0.146996 

-0.335611 

1.127617 

41 

0 . 666606 

-1.469497 

0.146996 

-0.285518 

1.105863 

42 

0.723848 

-1.456874 

0.146996 

-0.236273 

1.083745 

43 

0.777782 

-1.444335 

0.146996 

-0.187542 

1.060991 

130 


44  0.827561  -1.431955  0.146996  -0.138458  1.037087 

45  0.872381  -1.419472  0.146996  -0.088061  1.011463 

46  0.911495  -1.406547  0.146996  -0.035055  0.983411 

47  0.944210  -1.393024  0.146996  0.023027  0.951546 

43  0.969875  -1.37986 8  0.146996  0.091342  0.912385 

49  0.987862  -1.368515  0.146996  0.178655  0.361777 

50  0.997429  -1.365097  0.146996  0.303827  0.784388 


CD  =  0.030956  CL  =  0.713821  CM  =  -0.190635 


TRAILING  VORTICES  DATA 

M  X(M)  Y(M)  CIRC 

1  2.383780  0.232008  -0.000933 

2  2.333846  0.225080  -0.001625 

3  2.284404  0.216393  -0.002229 

4  2.235273  0.206695  -0.002784 

5  2.186347  0.196301  -0.003304 

6  2.137600  0.185376  -0.003797 

7  2.089004  0.174067  -0.004256 

8  2.040442  0.162550  -0.004687 

9  1.991957  0.150853  -0.005100 

10  1.943572  0.138989  -0.005481 

11  1.895155  0.127193  -0.005801 

12  1.346658  0.115584  -0.006100 

13  1.798126  0.104170  -0.006349 

14  1.749496  0.093078  -0.006559 

15  1.700793  0.082359  -0.006720 

16  1.651988  0.072087  -0.006839 

17  1.603077  0.062348  -0.006896 

13  1.554079  0.053196  -0.006916 

19  1.505019  0.044648  -0.006885 

20  1.455917  0.036758  -0.006795 

21  1.406791  0.029591  -0.006648 

22  1.357685  0.023172  -0.006443 

23  1.308651  0.017529  -0.006177 

24  1.259734  0.012684  -0.005852 

25  1.211021  0.008642  -0.005465 

26  1.162620  0.005401  -0.005014 

27  1.114706  0.002948  -0.004498 

28  1.067624  0.001210  -0.003917 

29  1.022530  0.000276  -0.003269 


TIME  STEP  TK  *  1.999990  TK  -  TKM1  =  0.200000 

ALPHA (T)  =  7.500000  OMEGA (T)  =  0.000000 

U(T)  =  0.000000  V(T)  =  0.000000 

NITR  VXW  VYW  WAKE  THETA  GAMMA 

0  0.939783  0.026143  0.188029  0.027811  0.156187E+00 

1  0.943367  0.030223  0.139770  0.031363  9.160749E+00 

2  0.948647  0.030031  0.139325  0.031699  j. 1607655+00 


CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =  2 

J  X(J)  Q(J)  GAMMA  CP(J)  V(J) 

1  0.997429  1.116135  0.160766  0.320772  -0.851401 

2  0.987862  1.110748  0.160766  0.221351  -0.907726 

3  0.969875  1.126034  0.160766  0.158555  -0.941336 

4  0.944210  1.150497  0.160766  0.116479  -0.962935 

5  0.911495  1.179177  0.160766  0.086724  -0.977640 

6  0.872381  1.210700  0.160766  0.065709  -0.987588 

7  0.827561  1.244766  0.160766  0.051593  -0.993866 


8  0.777784  1.281128  0.160766  0.043212  -0.997146 

9  0.723850  1.319423  0.160766  0.039694  -0.997912 

10  0.666609  1.359948  0.160766  0.040813  -0.996296 

11  0.606945  1.402723  0.160766  0.046325  -0.992422 

12  0.545769  1.448176  0.160766  0.056472  -0.986157 

13  0.484008  1.496484  0.160766  0.071440  -0.977363 

14  0.422591  1.547997  0.160766  0.091673  -0.965767 

15  0.362439  1.603142  0.160766  0.117358  -0.950894 

16  0.3C4449  1.662897  0.160766  0.151404  -0.931343 

17  0.249485  1.727721  0.160766  0.193674  -0.907603 

18  0.198363  1.798850  0.160766  0.247145  -0.376333 

19  0.151840  1.878118  0.160766  0.315652  -0.834967 

20  0.110606  1.967481  0.160766  0.404369  -0.778575 

21  0.075269  2.071099  0.160766  0.522025  -0.697329 

22  0.046353  2.194821  0.160766  0.679683  -0.571292 

23  0.024285  2.349161  0.160766  0.381013  -0.350465 

24  0.009388  2.539111  0.160766  0.937480  0.119060 

25  0.001884  2.420488  0.160766  -0.772358  1.331331 

26  0.001884  -0.236994  0.160766  -4.940613  2.437122 

27  0.009387  -1.441366  0.160766  -3.294956  2.071290 

28  0.024284  -1.678374  0.160766  -2.145366  1.771569 

29  0.046353  -1.735029  0.160766  -1.575494  1.601988 

30  0.075269  -1.740161  0.160766  -1.248166  1.495590 

31  0.110606  -1.726694  0.160766  -1.035354  1.421363 

32  0.151840  -1.705142  0.160766  -0.384674  1.367020 

33  0.198363  -1.680851  0.160766  -0.770207  1.323621 

34  0.249485  -1.655530  0.160766  -0.678837  1.287744 

35  0.304448  -1.630341  0.160766  -0.602871  1.256973 

36  0.362436  -1.606176  0.160766  -0.537036  1.229567 

37  0.422588  -1.582802  0.160766  -0.478581  1.204600 

38  0.484005  -1.560513  0.160766  -0.425234  1.131273 

39  0.545766  -1.539373  0.160766  -0.375344  1.158993 

40  0.606941  -1.519384  0.160766  -0.327563  1.137220 

41  0.666606  -1.500277  0.160766  -0.231143  1.115663 

42  0.723848  -1.482141  0.160766  -0.235015  1.093871 

43  0.777782  -1.464664  0.160766  -0.188536  1.071527 

44  0.827561  -1.448072  0.160766  -0.140512  1.048044 

45  0.872381  -1.432248  0.160766  -0.089803  1.022803 

46  0.911495  -1.417014  0.160766  -0.035089  0.995022 

47  0.944210  -1.402392  0.160766  0.026050  0.963233 

48  0.969875  -1.389331  0.160766  0.098463  0.924430 

49  0.987862  -1.379196  0.160766  0.190332  0.873000 

50  0.997429  -1.379105  0.160766  0.319847  0.795184 

CD  =  0.022037  CL  =  0.709635  CM  ■  -0.185327 


TRAILING  VORTICES  DATA 

M  X(M)  Y(M)  CIRC 

1  2.923984  0.318366  -0.000933 

2  2.873369  0.311788  -0.001625 

3  2.823534  0.302817  -0.002229 

4  2.774192  0.292439  -0.002784 

5  2.725182  0.281064  -0.003304 

6  2.576492  0.268880  -0.003797 

7  2.528086  0.256054  -0.004256 

3  2.579755  0.242831  -0.004637 

9  2.531597  0.229197  -0.005100 

10  2.483748  0.215114  -0.005481 

11  2.435886  0.200928  -0.005801 

12  2.387918  0.186791  -0.006100 

13  2.339993  0.172653  -0.006349 

14  2.291922  0.158710  -0.006559 

15  2.243772  0.144992  -0.006720 

16  2.195477  0.131596  -0.006839 

17  2.146968  0.118650  -0.006896 

18  2.098253  0.106217  -0.006916 

19  2.049394  0.094288  -0.006885 

20  2.000362  0.082951  -0.006795 


21  1.951119  0.072314  -0.006643 

22  1.901708  0.062412  -0.006443 

23  1.352145  0.053298  -0.006177 

24  1.302427  0.045045  -0.005352 

25  1.752599  0.037677  -0.005465 

26  1.702674  0.031252  -0.005014 

27  1.652678  0.025824  -0.004498 

28  1.602603  0.021481  -0.003917 

29  1.552431  0.013390  -0.003269 

30  1.501920  0.017193  -0.002553 

31  1.450583  0.016109  -0.002761 

32  1.378677  0.012411  -0.005504 

33  1.253439  0.007139  -0.007734 

34  1.094864  0.003008  -0.009244 
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